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Abstract 


The dynamic response of (i) a simple beam, and (ii) a 
Single bay portal frame supporting an unbalanced rotor is 
investigated in detail using a computer simulation. The 
method of solution for transient response is paced on direct 
(step-by-step) integration of the system equations of motion 
through application of a finite time element reccurence 
scheme. Finite, Timoshenko beam elements are used ih 
modelling both mass and elasticity as distributed parameters 


of the supporting structures. 


Two modelc scot the “forcing function, due to rotor 
unbalance are considered, assuming a rigid rotor. shaft 
Supported on: (i) rigid bearings; and (11)-011-film journal 
bearings. The frequency of the forcing function is time 
dependent to simulate transients at start-up or shut-down 
operations. The discussed Supporting structures are 
low-tuned relative to the EOCOL operating speed. 
Consequently, the time-dependent frequency of the excitation 
force passes through at least one critical frequency of the 


foundation system. 


The results of the transient analysis are presented in 
the graphical form and discussed. The analysis 1s based on 
non-dimensional parametric studies of the system. The study 
shows that the maximum amplitude of vibration of low-tuned 
Supporting structures is highly dependent on rotor 


acceleration rate through the critical (i.e. natural) 
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frequency of the foundation system. When compared to the 
results of classical steady-state analysis, the maximum 
amplitude obtained through the transient analysis is greater 
in magnitude and itS position is shifted. These results 
indicate, that in the case of low-tuned structures, the 
transient analysis should be considered as a_e standard 


procedure in the present-day design practices. 


The study suggests that the numerical procedure, and 
the computer program developed for this purpose, are useful 
for the transient analysis of the models of foundation- 
structure interacting with rotating machinery. The method of 
solution is general for any model of forcing function, and 
the computer program can easily be extended to handle more 


complex, two or three dimensional structures. 
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5.36 Effect of rotor unbalance on system dynamic 
response for rotor speed parameter a=0.8, for 
model with (1) journal bearing, and (2) rigid 
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5.37 Oil-film force magnification factor versus time, 


for two different values of rotor unbalance 
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1. INTRODUCTION 


1.1 Background Information 


It is well known that the life, efficiency and overall 
performance of all kinds of rotating machinery depend 
largely on the vibration of the machine and its supporting 
Structure. An extensive investigation in the field of 
machinery vibration has been carried out for many years with 
numerous results of experimental and analytical works being 
published. Remarkable progress has been achieved in this 
field, especially during the past twenty-five years with the 
development of powerful, fast computers and new numerical 


techniques. 


The demanding deSign requirements placed on modern 
rotating machinery have introduced a trend towards building 
units of larger size, higher speed, increased power and 
efficiency. As a result the task of minimizing and 
controlling the system vibration levels has become even more 
important. The trend to build larger rotary machines which 
rest on masSive reinforced concrete or steel pedestals has 
brought changes into traditional design of machine 
foundation system and also created many problems in 
Structural dynamics. Some of those problems are closely 
associated with the phenomenon of dynamic interaction 
between an elastic structure and a deformable machinery. 


While extensive research has been carried out on the rotor 
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dynamics, relatively little work is done on the 
machine-foundation interaction. It should be mentioned, 
however, that this problem has been recently attracting 
considerable attention from both the industrial and research 
communities. The dynamic response of the machine- foundation 
interacting system, in general, depends on design criteria 
and characteristics of the specific class of machinery and 


its supporting structure. 


Figure 1.1' shows a schematic diagram of a large 
PEO Generator set, sand F1qg,eni.2,1 sits »foundatren sblock. 
(This foundation block was resently built. Its design and 
innovative features are outlined by Ellyin [1]?). There are 
four distinct elements, namely: (a) rotor, (b) bearings, 
(c) bearing pedestals, and (d) foundation block. They act as 
one system, together responding to rotor acceleration, 
critical speed, external excitation forces, et cetera. It is 
obvious that any change in the subsystem will affect the 
vibration levels of the other parts and, consequently, the 
dynamic response of the system as a whole. Therefore, 
a basic model for the dynamic analysis of the machine- 
foundation interacting system should include all those 


distinct elements. 


The rotor, bearings and bearing pedestals are designed 


by a turbine manufacturer. The foundation block is usually 
‘First number denotes Chapter. Figures are grouped at the 
end of each chapter in the order they are referenced to. 
*numbers in square brackets designate Reference at the end 
of the thesis 
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designed by another organization using data and requirements 
provided by a machine manufacturer. The foundation block is 
a very important part of the system. It provides not only 
mass, flexibility and damping but also the essential 
coupling between other parts. There iS a trend to _ build 
flexible foundations which are "low-tuned" (or 
"under-tuned") relative to the machine speed. In Such cases, 
the machine service speed is at least higher than the first 
natural frequency of the Supporting structure. Consequently, 
as the rotor is brought up in speed during start-up 
Operation, the frequency of the forcing function? (i.e., 
exciting force) passes through one or more critical (i.e., 
resonant) frequencies of the foundation system. This can be 
best illustrated by a resonance diagram for an idealized 


Single-degree of freedom model, as shown in Fig. 1.3. 


Due to machine-foundation interaction, a turbomachinery 
designer must analyze the dynamic performance of the 
integrated system. The problem is extremely complex, and 
even with analytical and numerical tools available today, it 
constitutes a formidable task requiring a certain degree of 


model simplification. 


In the present study the problem is approached from the 
foundation system designer's point of view. Accordingly, the 


analysis focuses on the dynamic response of a supporting 


*The sources of the exciting force are numerous. Here only 
the most common one (that due to rotor unbalance) is 
considered 
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structure excited by a force transmitted from the rotating 
machinery. Consequently, the difficulty now lies in defining 
aeemodel® of Sthe® forcing? function®@@which approximates the 


exciting force with reasonable accuracy. 


Massive foundation blocks supporting modern turbmachin- 
ery are generally constructed from simple beam and frame 
elements (Fig. 1.2), constituting the main load transfer 
members of the supporting structure. To appreciate the 
dynamic behavior of the whole system, it is of great 
importance to (fully) understand the dynamic response of 
these simple elements. While there is an abundance of 
literature on steady-state modai analysis of beam and frame 
models, in contrast very few works are published on 


transient analysis of these elements. 


1.2 Scope of Study 


The dynamic response of a simply supported beam 
Subjected to a force of time-dependent frequency has been 
invesStigated+by’=Suzuki" [2-4]vand* Victor ©& Ellyin:s [5]. An 
analytical approach to solving model equations of motion 
enabled the authors to carry out thorough parametric studies 
and to draw important conclusions. However, the proposed 
methods of solution limit their practical application to 
only: (a) simple structures, and (b) certain models of 
exciting force, which could be expressed by relatively 


Simple analytical functions. 
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The present study attempts to fill a gap by basing 
method of solution on ae finite element technique, thus 
making it suitable for more general models. 

The main objectives of the present work are formulated as 


follows: 


1. To develop a method of solution and a computer 
program, based on the principles of finite element 
technique, suitable for transient vibration analysis 
of different supporting structures subjected to any 


type Of gforcing Ehunchuon: 


2.) 1O verify the “method of ° solution and test the 
computer program, by applying it to transient 


analysis of a model, already investigated by others. 


3. To perform a transient analysis of a portal frame, 
as a model of foundation structure, supporting an 


unbalanced rotor mounted on rigid bearings. 


4, To develop a model of a forcing function, assuming 
anncunbalkancede srotor esmounted)] ron] oilsirlmajournal 
bearings, and analyze dynamic behavior of a_ simple 


beam subjected to this excitation. 


To achieve these objectives the study is divided into 


four stages summarized below: 


Phase 1. Firstly, a general n-degree of freedom 


structure and its system of second order dynamic equations 
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of motion is considered. Next, a finite element technique 
(applied to the time domain) is used to derive a finite-time 
formulation. Then, an initial value-transient problem is 
assumed and resulting reccurence scheme derived (Chapter 2). 
Finally, a computer program is developed which employs the 
reccurence formulae to advance the solution step-by-step in 


time (Chapter 4 and Appendix A-3). 


Phase 2. A model (Model 1 -"Beam Supporting an 
Unbalanced Rotor Mounted on Rigid Bearings") is defined 
USecElonioec)..thesrorcing functions proposed “by = Victor | & 
Ellyin [5] is adopted (so that readily available results can 
be used) to verify the method of solution and computer 
program. Results of detailed parametric studies of the model 


are presented and discussed (Section 5.1). 


Phase 3. A model (Model 2 -"Single Bay Portal Frame 
Supporting an Unbalanced Rotor Mounted on Rigid Bearings") 
is described (Section 3.3). Two components (in horizontal 
and vertical directions) of thesexciting force, idue to rotor 
unbalance, are considered. Results of numerical analysis, 


discussion and conclusions are presented (Section 5.2). 


Phase 4. A model (Model 3 -"Beam Supporting an 
Unbalanced Rotor Mounted on Oil-Film Bearings") is developed 
(Section 3.4). Short bearing approximation, as suggested by 
Ocvirk [6], and non-linear performance of journal bearings, 
Holmes 7 Sel are assumed. Results, discussion and 


concluding remarks are presented (Section 5.3). 
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A final chapter (Chapter 6) links the separate findings 
together and provides an overall conclusions and 


recommendations resulting from the study. 
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Fig. 1.1 Model of rotor and foundation system 
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Fig. 1.2 Sketch of turbo-generator foundation block 
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Fig. 1.3 Single-degree of freedom rotor-foundation model and 
its resonance diagram 
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2. FINITE TIME FORMULATION TO INITIAL-VALUE TRANSIENT 


PROBLEM 


2.1 Preliminaries 


The transient response analysis of a structural system 
generally involves discrete modelling of the structure, 
whereafter the continous physical structure iS approximated 
by a discrete n-degree of freedom model. Regardless of the 
Spatial discretization scheme employed, the resulting set of 


System dynamic equations of motion (in matrix form) becomes: 


[Mliq}+([Cliq}+[Kliq}={F} (20.1) 


where‘: [M]®, [C] and [K] are the (nxn) mass, damping and 
stiffness matrices; {F} is the (nx7) external load vector, 
ands 1qt,-10}" and wyg) are the (nx?) nodal” “acceleration, 


velocity and displacement vectors, respectively. 


The set of ordinary differential equations (2.1) can be 
integrated forward in time to generate a transient solution. 
There are a great variety of direct time integration 
schemes, which fall into one of two major categories (or 
their combination), namely: explicit procedures and implicit 
procedures. Both methods have substantial advantages for 
certain classes of problems and disadvantages for others. 
Each specific scheme has different accuracy, stability 


*A complete list of notation used is given in Appendix A-1. 
5In some expressions it is appropriate to simplify matrix 
and vector notation. Throughout the text two notations are 
used interchangeably, e.g. [M]J=M, {q}=q. 
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characteristics and numerical efficiency. (An excellent 
review of the most widely used schemes and their suitability 
for various engineering problems is given by Donea, et al. 


aR 


The time integration scheme developed for the purpose 
of the present analysis is based on a finite element (in 
time) approximation technique. A_ finite time element 
consists simply of a fixed time interval, At, which can be 
treated as a standard, one-dimensional finite element with 


Ewounodeseat t=t, and t=t> (At=r=t.-t,). 


The general finite-element discretization process can, 
therefore, be applied to the time domain. First, the 
interval At is treated as a finite domain of time, anda 
finite time formulation is derived for the full original 
matrix equations (2.1). This formulation relates the values 
of vectors {q}, {q} and {q} at time t, and t2. As the time 
dimension is of an infinite extent, the solution for the 
initial value-transient problem is obtained by repeating the 
calculations for subsequent finite domains of time with new 
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2.2 Derivation of Reccurence Formulae 


Consider a two-node element with three degrees of 


freedom per node: 


q1 ve 
gree == eee a2 
q: t; 7 t2 q2 


This gives six "time-wise” degrees of freedom per element. 
The problem now is to find a function q(t) which satisfies 
the equations of motion (2.1) and the following boundary 
conditions: 


q(t,)=q, a(t,)=a; qa(t,;)=a, 
(252) 


q(t2)=q2 q(t2)=q:2 a(t2)=q; 
In order to approximate the function q(t), satisfying the 
conditions (2.2), at least a fifth-order polynomial in time 
is required, as follows’: 

: Ge) = GQ, 7a2t+¢s0+34a,C¢tast-ta6t? 
(223) 

=e ent iC (eect ee cm cee te. a 
or rather, remembering that q(t) represents a vector (nx?) 
of nodal displacements with each component approximated by 
the same polynomial, a set of such polynomials is required, 


as follows: 


fleece lh. stele tel. ter tat 


= [@]" {a} 


‘The following formulation was provided to the author by the 
thesis supervisor, Dr F. Ellyin. 
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where: 
I - unit matrix of order (nxn) 
[@] - matrix of time terms (6nx6) 


{a} - vector of undetermined coefficients (6nx7) 


The first» and second derivatives of q(t) are obtained from 


eq. (2.4) as: 


[]" {a} 
Be ab225) 
[$] {a} 


qe) fOe*E 2EPr* '3t22* 42 re *5t4rifa} 


q(t) L0e%0 21 6tie tet he 20t hia 


The (6nx1) coefficients {a} can be evaluated in terms of 
nodal variables by matching the displacements, velocities 


and accelerations at the nodal points t,=0 and t2=r, which 


leads to; 
{qi} I 0O 0 0 0 0 varst 
{qi} 0 I 0 0 0 0 {a2} 
{qi} O10 6-24 0 0 0 {a3} 
{q,it= = (256) 
{q2} oer? 7 tt Fal Tot ren fa,} 
{q2} 0 Fpe2xliearai mane E5rct {as} 
{q2} 0 0 21 6rI 12721 20722 {ag} 
(6Nx17) (6nx6n) (6nx7) 
or, 
{qn t=[T]{a} (237) 


Therefore, 
faj=l To ice (28) 


a 


i te Shae %, alley 
gy 8}. 499 ae ic 


| a io, 0 
eangiatPieon (mel ef 


soi teietes 


wheres 
27,27 0 0 0 0 0 
0 2771 0 0 0 0 
0 0 TAY 0 0 0 
ewe i / 27° 


-20r°®I -127r’7I -3r*I 20r§SI -8r’I r®I 
S077) Vl6rol, 3771 =307-1 147°! -27r7 1 
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Hence, substituting eq. (2.8) into eqs. (2.4) and) (2.5) 
giveS an approximation for displacements, velocities and 


accelerations in terms of nodal variables {q,}, as follows: 


q(t) = [@][T]-'{q,} = [N]{qn} 
q(t) =~ [6] '[T]-'{q,} = [N]{qn} (220) 
Q(t) =~ [6] [T]-‘{q,} = [N]{qn} 


where: 


[N]=(]' [T]-'=[N, Nz N3 No Ns Neo] 


is the matrix of interpolation or shape functions, which in 
this case are fifth-order Hermitian polynomials [10]. 
Substituting eq. (2.9) into the equations of motion (2.1) 


gives 
[MJ(N]{q,}+(C](INJ ig, }+(KI IN] {q,}-{F}#0={R} (2750) 


Note that upon substitution of the approximated functions 
(2.9), the equations of motion will generally not be 


Satisfied. In fact, there is some residual {R} left. 
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Now, the Galerkin (weighted residual) method can be 
employed to find / unknowns of the vector {q,} in such a way 
ehatethe residual, “tRipb over ‘the! “entire domain; +r, is 
minimized with respect to the shape functions, N;. This 
process can be expressed in the form of an integral equation 


as 


T 
JN, iR}dt = 0 
Ory, 


SN, (IMIIN] (qn }+{C] (NI fa, }+(K]IN] {q,}-{F})dt=0 (etn) 


Equation (2.11) represents j simultaneous equations in the j 
unknowns. In the case of an initial value problem, q,, q; 
and q; are assumed known and the above equation serves to 
approximately determine q2, q2z and qz2. Quite generally, 


however, equation (2.11) can be rewriten in the form 
SIN] (IMI(N] {qn 3+(C](N] fa, }+(K][N] {q, }-{F})dt=0 C2712) 
Substituting relationships (2.9) and [N]'=[T-']’ [@] into 
(2.12) gives 

[v-}" s" ({L@] (MI ($]' +Le](C](o]"+(S](K][8]" FIT] fq,} 


-[@J{F})dt = 0 C2ii3)) 


Carrying out simple integration and matrix multiplication 


leads to 


Zi tbc ele aoe ate) (2.14) 
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where: 
I {P,} 
tI {P.} 
Talc {P3} 
Poe ee iE vata ies (22-45) 
cet {P,} 
jw | {P,} 
or, 
T 
{P,} = S deemed Ol EB Keke P= 6 (22516) 
0 -M 0 0 M 0 
M 0 0 -M ™ 0 
™ ir?M soT?M -7M 277M ) +¢7°?M 
CZ e= G2eali7s) 
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Introducing 


Cee a oral, 


into equation (2.14), the general relationship can be 


written as 
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: {q:} {P,} 
Pees: Mees ie a) 
{qi} fea 
coerce ne ec cece ec eoe | [ooee | = | ovee G22. 1) 
{q2} {P,} 
P7eeelee CZ. no taet (oan 
: {qa} {P,} 
(6nx6n) (6nx1)  (6nx7) 


For an initial value problem, only part of the general time 
formulation (2.21) is required to determine q2, q2z and q2. 


Thus, 


{qi} {q2} {P,} 
244 toa + Z12 {qo}|} = |{P2} C222) 
{qi} {q2} {P3} 


Applying simple matrix algebra and rearranging terms in eq. 


(2.22) leads to the formulae; 


{q2} ae {P,} {q,} 
iq2}|] = Z42 | Bois) bare 244 faa} (2.23) 
{q2} {P3} {q:} 


In the foregoing it has been assumed that the domain of 


the approximation corresponds to interval At. Equation 
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2952'3;) represents a reccurence relation between two 
Sucessive values of vectors q(ttAt), q(ttAt),q(ttAt) and 
q(t), q(t), q(t). This relation will be used in a computer 
program for a sequence of such domains in a_e step-by-step 


manner to advance a solution in time. 


The step-by-step reccurance scheme just derived is 
conditionally Stable. This means that for numerical 
Stability (convergence of solution), the time step At cannot 
exceed some critical value. This critical value, in general, 
is "controlled" by the system's highest natural frequency 
(wm) and the requirement for stability is At<a/w,. A value 
of constant "a" depends on stability characteristics of a 
Specific scheme used. Stability analysis for several widely 
used schemes and resulting analytical expressions for a's 
are given in reference [11]. Such analysis for the present 
scheme (due to the high order interpolation used) proved to 
Demea friculy.. site ShOULG  SILfice to “say that. from the 


practical experience of this work a=1.4. 


In general, the dynamic equations of motion (2.1) can 
be non-linear if, for instance, the damping matrix is 
frequency dependent or, as in the case of large deformation, 
the stiffness matrix depends on displacements. However, in 
the present analysis, only linear systems are considered. It 
is, therefore, assumed that the system [M], [C] and [K] 
matrices are not time dependent and consequently neither are 
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non-linear system, a small amount of computing time and 
Simple program organization is required. Matrix [2,2] has to 
be inverted only once and then, for each time step, only 
integration (to evaluate components P,, eq. 2.16) and simple 


matrix multiplication is involved. 


3. PHYSICAL MODELS DESCRIPTION 


3.1 Preliminaries 


In the present analysis the dynamic response of two 
Simple models of support structure subjected to the action 


of an unbalanced rotor are examined, namely: 


1. a beam simply supported at its ends 


2. a Single bay portal frame with both legs clamped 


Distributed inertial and elastic properties of these 
structures are modelled, using the consistent mass and 
stiffness matrices of the finite element approximation [11]. 
A number of finite beam elements have been proposed in the 
literature [12-21]. They are based on_- either: (i) 


Bernoulli-Euler or (ii) Timoshenko beam bending theories. 


The classical Bernoulli-Euler’ theory of flexural 
vibration of beams considers only the lateral inertial and 
elastic forces due to bending deflections. This theory is 
Satisfactory for dynamic analysis of long, thin beam models 
with aesmall slenderness ratio, k/G (i.e. VI/A/L)< However; 
large rotating machinery support systems are usually 
composed of short stubby beams. In such cases, the secondary 
effects of the shear deformation and rotary TAerETs of the 
cross-section of the beam cannot be neglected. These 
secondary effects are included in the Timoshenko beam 


theory. Therefore, the finite elements based on this theory, 
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known as the Timoshenko beam finite elements, are used in 


this study. 


Two different models of the Timoshenko beam finite 
element are used in the present work, namely those: (a) 
proposed by Davis, et al [13], and (b) by Akella & Craggs 
[21]. The nodal variable u (axial displacements) is added to 
the original beam elements to enable axial-flexural coupling 
in the framework models (see Fig. 3.1). For easy reference, 
mass and stiffness matrices for both these elements are 
given in\explicit form in Appendix A-2. The generalized 
coordinates at each node of element (a) represent: w- 
transverse displacement, @- cross-section rotation, u- axial 
displacement, and for element Gent wr transverse 
displacement, ¢- bending slope, AyW- shear slope, hence shear 


force, I¢'- bending moment and u- axial displacement. 
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Fig. 3.1 Timoshenko beam elements 
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The aforementioned element matrices are used in the 


structure discretization process as follows: 


1; >cA°*modelcof thesconcerned structure is divided into a 
number of elements, for which each of the element 


mass [EM] and stiffness [EK] matrices are evaluated. 


2. A standard finite element assemblage process is 
employed to form global mass [GM] and stiffness [GK] 


matrices. 


3. Global matrices are modified to incorporate system 


constraints and boundary conditions. 


Upon completion of this proces, a continous structure is 
approximated by using the spatially discrete m-degree of 
freedom model with (nxn) system mass [M] and stiffness [K] 


matrices. 


Borormostenpracticala stnuctures; ethetrexact form of 
damping is unknown, and damping properties are frequency 
dependent. However, in order to take advantage of the 
explicit time stepping scheme, it is necessary to evaluate 
the damping matrix explicitly. One procedure for defining a 
system damping matrix is to employ a particular form of 


proportional damping, given as 
[Cl=alM)e seb ik) G37 an) 


where, @ and b are proportionality constants. Such damping, 


known as "Rayleigh damping" [22], is used in the present 
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analysis. 


As mentioned in the Introduction, a self-excited 
vibration of rotating machinery could be generated by many 
different sources of excitation (e.g. rotor imbalance, 
internal Eri ction. misalignment, gyroscopic moments, 
hydrodynamic fluid-film forces, et cetera [23]). In the 
present study, only the most common source of excitation, 


i.e., the rotor unbalance, is considered. 


A force transmitted from an unbalanced rotor to its 
foundation iS an extremely complex function of many 
parameters. In general, it depends on rotor geometry and 
PhexebDiity, and on mounting details, that is 
characteristics of bearings, bearing seals, bearing 


pedestals, and so on. 


In the present work the following simplifying 


assumptions (general for all models) are made: 


1. Symmetric rotors are represented by a lumped mass 
(2Mr ) located centrally between two identical 


bearings. 


2. Rotor shafts are massless and absolutely rigid. 


3. Bearing pedestals are neglected. 


4, Loads are only in one plane. 
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25 
3.2 Model 1 - Beam Supporting an Unbalanced Rotor Mounted on 
Rigid Bearings 
3.2.1 Forcing function 


To further simplify the study, rigid bearings are 


assumed for the model shown in Fig. 3.2. 


Fig. 3.2 Beam supporting an unbalanced rotor 


Taking all the assumptions from the previous Section 
into consideration, the forcing function for this model can 


be expressed by the following equation: 
F,(t)=-Mr-d?[e-sinQ®(t)thty,(t)]/dt? OS S23) 


Carrying out the deriviation leads to 
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FY (t)= Mree(Q?sinQ-YcosQ) - Mreym ig ey 
or 
F,.(t)=F,(t)-Mr«y-a(t) (3.4) 
where: 
By - a vertical component of an inertia force due to 


rotational motion of an unbalanced rotor 
Mr«ym - an inertia force due to rotor acceleration in 


its linear motion 


As in reference [5], it is assumed that the torque applied 
to the rotor by the power source is constant during the 
acceleration time (T,). This means that in an ideal 
Situation (with no vibration) the rotor angular acceleration 
would also be constant. In reality, however, part of the 
torque is absorbed by the vibration. Because of the build up 
of vibration, the torque associated with the unbalance 
reduces the rotor angular acceleration. Thus, the angular 
acceleration %, velocity ®, and travel 2, are assumed _ to 


vary with time according to the formula: 
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27 
3.2.2 Finite element model of supporting beam 


A typical transversal beam of the reinforced-concrete 
foundation block (shown in Fig. 1.2) is chosen for the 
dynamic analysis. Its geometric parameters and material 


constants are given below 


L= 300.0: in v= 0.15 K= 0.85 
A= TAS0 20° in" p= (4629 p/EtC? 


I= 9266400.0 in” E= 5500000.0 lb/in? 


Fig. 3.3 A typical transversal beam of the foundation block 


The high order finite element (TM544) is used to 
approximate inertial and elastic properties of the continous 
beam. This finite element is chosen for’ the numerical 
analysis because, in comparison with other available 
Timoshenko beam elements, it iS more accurate and converges 
faster per degree of freedom for beams with k/L>0.05. The 
beam selected for the analysis (Fig. 3.3) has a k/L ratio 


equal to 05712. 
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The accuracy of a finite element approximation depends 
on the number of elements used in the discretization 
process. In general, the greater this number, the more 
degrees of freedom considered and the better accuracy. 
However, at the same time, the size of the problem and 


consequently the cost of numerical analysis is greater. 


In order to assess a minimal number of finite elements 
required for beam approximation (sufficient for the purpose 
of this project), a preliminary study was carried out. 
Natural frequencies and dynamic response of the beam 
idealized by 2 elements and by 4 elements (which gives 9 and 
19 D.O.F of the constrained systems respectively) were 
calculated and compared. No appreciable difference in 
dynamic response but dramatic increase in the cost of 
solution was observed in the second case. This can be 
explained by the fact that the response is dominated by the 
frequencies of lower modes which are predicted with 
sufficient accuracy by the 2-element model. On the other 
hand, the dimensions of the finite element matrices, which 
determine the size of the problem, increased from (27x27) in 
the first case to (57x57) in the second. Moreover, since the 
time step required for numerical stability is controlled by 
the highest mode frequency, for the 4-element model this 
time step has to be reduced considerably. Taking the above 
into consideration, the 2-element model approximation is 


used for further analysis of the beam. 
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3.3 Model 2 - Portal Frame Supporting an Unbalanced Rotor 


Mounted on Rigid Bearings 


Seas FOrcing function 


All the simplifying assumptions discussed in the 
previous section are also applied to Model 2, of which a 


schematic diagram is shown in Fig. 3.4. 


Fig. 3.4 Frame supporting an unbalanced rotor 


In this model, two components of the forcing function, 
are considered, in horizontal and vertical directions. 


Applying the same reasoning aS in Model 1 leads to _ the 
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expressions 


F,(t)= Mr-e(Q?cosQ+QsinQ) - MreXm 
a6) 
Fy(t)= Mree(Q?sinQ-QcosQ) - Mreym 


where 2, 2 and &% are assumed to vary with time according to 


20S. 03.5) 


3.3.2 Finite element model of supporting frame 


A typical frame of the turbo-generator foundation block 
(shown in Fig. 1.2) is chosen for numerical analysis as a 
Supporting structure in Model 2. Geometric parameters of the 
frame columns are given in Fig. 3.5. Material constants and 


beam geometry are identical to those in Model 1 (Fig. 3.3). 


Column geometric parameters: 


L= 600.0 in A= 3980.0 in? T= 895500.0 in* 


Fig 3.5 A typical frame of the foundation block 
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A preliminary study done on this model showed no 
advantage in using complex Timoshenko beam elements for 
idealizing the frame. In fact, the use of simple elements 
manifests faster convergence per degree of freedom. This is 
so because the columns of the frame selected for numerical 
analysis are fairly slender (k/L=0.025), while the high 
order elements prove to be superior for shorter stubby beams 
in which shear and rotary inertia are particularly 
important. The simplest Timoshenko beam element (with 
centroidal displacement and cross-section rotation as nodal 
variables) is the most suitable for the analysis of complex 
Structure and the most commonly used in general purpose 
computer packages [20-21]. As a result of the preliminary 
study, the 6-element model approximation was chosen for the 
dynamic analysis of the frame, giving 15 D.O.F. for the 


constrained structure. 


In the present finite element model of the frame, there 
is geometric misrepresentation of the structure at the 
corners where two beams are joined at right angles (see Fig. 
3.5). However, this misrepresentation can be avoided in 
general analyses of frameworks by developing beam finite 


elements with special ends [13]. 
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3.4 Model 3 - Beam Supporting an Unbalanced Rotor Mounted on 


Oil-Film Journal Bearings 


3.4.1 Introduction 


Hydrodynamic bearings play a very important part in 
rotating machinery vibration problems. Their three major 
functions are: (a) to support loads (static and dynamic), 
(bopitos control@rotor position, and (c) to provide stiffness 
and damping. Bearing film damping has a dominant effect on 
rotor vibration. The importance of correct selection of the 
bearing and lubricant parameters, depending on design load 


and speed, cannot be exaggerated. 


The dynamics of rotor-bearing systems has become of 
great importance due to the technological advancement in 
modern rotating equipment, and has been studied extensively 
6047.0 24--34|, The field is very broad and complex. 
Researchers uSually devise simplified models and focus their 
studies on particular problem such as, system instability, 
response to shock or unbalance excitation, effect of rotor 
and support flexibility and/or damping on system response, 
et cetera. The vast majority of papers reviewed by the 
author, in relation to this project, consider the system 
behavior about its equilibrium position, that iS assuming 
constant rotor speed. In fact, the complete transient 
analysis of rotor-bearing system during the whole period of 


Start-up until the steady state operation does not seem to 
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have been reported in the literature. 


In the following section, a simplified analysis is 
presented, which leads to determination of the dynamic 
Oil-film forces induced in the journal bearing by an 
unbalanced accelerating rotor. The resultant of these forces 
(in vertical direction) is then considered in the original 
problem aS a component of forcing function exciting rotor 


Supporting beam. 


The analysis is based on the following simplifying 


assumptions: 


1. An ideal, symmetric rotor allows to consider’ the 


motion of one journal only. 


2. The Reynolds' equation (lubrication equation) for an 


incompressible lubricant is applicable. 


Se sViSCcOSity 1S constant throughout; the, poibsirim, and 


does not change with time. 


4, Oil-film pressure at the ends of the bearing is 


atmospheric. 


5. Flow of lubricant in the bearing axial direction can 


be neglected. 


6. Only positive pressure contributes to the bearing 


dynamic film forces. 
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3.4.2 Dynamic analysis of journal bearing 


Consider a journal-bearing system shown schematically 


Oj WA a Se PRS Se) 


C - radial clearance 5%, 
e = e,/C - journal centre eccentricity 
h = C(1 + ecos@) - oil film thickness 


Fig. 3.6 Journal-bearing system 


Consider first the case when shaft angular velocity w 
is constant. The journal is in motion due to the action of 
four external forces, namely: static load (weight) WwW, 
centrifugal force due to rotor unbalance P,, and 
hydrodynamic oil-film forces P{ and P%. In analyzing the 
motion of the journal, it 1S convenient to employ polar 


coordinates rather than Cartesian. (The polar coordinates 
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used are, the eccentricity ratio e and the attitude angle 
w). Assuming that the hydrodynamic force P% acts in the 
Positive ..e-direction si se. opposite ytosthat, ,shown pin (Fig. 
3.6), the equations of the journal centre motion may be 
written 

PS +WcoswtP,cos(wt-wy)=MrC(e-7e) 

P$-Wsinwt+P,sin(wt-wW)=MrC(ewt+2ey) a 
The set of equations (3.7) can be easily solved using one of 
many available numerical schemes. The only problem is, that 
it requires evaluation of the dynamic film forces (Pi, P$) 
at each time step. These forces are determined by 
integrating oil-film pressure distribution around the 
journal circumference. The pressure distribution, in turn, 
is obtained by integrating Reynolds' lubrication equation 
which, in general, has to be done numerically. Considering 
the fact that the time step required for solution of eqs. 
(3.7) has to be small enough to give 100 increments per 
evclesoresshnattawmotion | [32]) > sthe wcalculation; vot “bearing 
forces for relatively long transient period can prove to be 


very expensive in terms of computer time. 


The Reynolds' equation for a plain cylindrical bearing, 
whose journal centre position is time dependent, is usually 
Written in) the form [35]: 
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Holmes [7], reports that his experimental works support the 
use of short bearing approximation (due to Ocvirk [6]) in 
describing the dynamic behavior of modern turbine bearings. 
Ocvirk showed that for short bearing (L/D<1) dp/d0<<dp/dz, 
which allows to neglect the first term in equation (3.8). 
Thus 


d(h?0p/dz) /az=6u(w-2y) dh/d6+12uCecosé@ (3.9) 


This greately simplifies the problem, since the reduced 
Reynolds' equation (3.9) can be integrated analytically. 
Integrating twice with respect to z and applying boundary 


conditions: p=0’ at z=0 and at z=L yields 
D=3uz\z-L) P2ecosé=(w-2y) esind][C7(1+ecosé@ ) J= (3.10) 


The oil-film forces, P% and P$, can now be obtained by 
integration as follows: 
L 9, 
Pi = Rf J, pcosédédz 
8, 
L 9, Carr) 
P$ = RJ J, psinédédz 
0 “8 
Integration of eq. (3.11) is Straight forward provided that 
the limits of integration 9, and 962 are known. In his 
earlier work, Holmes [7], suggests that these limits lie at 
the points were the pressure, p, becomes zero (see Fig. 3.6) 


and can be determined from eq. (3.10) as: 
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According to his later paper [8], much better agreement with 


the experimental results is achieved by putting 
6kn=w0 cand &keswme 1 G3.=1.30) 


CanuyingesOut integration of eqs. (3,14), with the limits of 
integration given by (3.13) yields 
Pau oak £2C faskn Git Beis) eth? (a=2y be 401 Se 4) PoSGi—es) ~ 247 
(3.14) 
P2= § (eb RA2C A phYede tmei(o=2yi) ( -a7:) B32] (1 Fe" 7 2? 
Under static conditions (i.e. in an ideal case of perfectly 
balanced rotor) the journal centre takes an equilibrium 
position (€o9,W.o), which depends on bearing geometry (L,R and 
oC), iMbwicantmaviscosityee(n)iperstatic load (Ww) andtangular 
shaft velocity (w). This position is given by the well known 


relationships [35]: 


GG ea) tleptt6estr= (mes ri Aye! 63-115)) 


Wostan tha ie 3) at /4eg] (3.16) 


where S,,, known as modified Sommerfeld Number, is defined 


as: 
Sm=uHL?Rw/4C7W Can aie) 
Relationship (3.17) can be rewritten 
(uL?R/2C? )=2S ,W/w C3382) 


It is apparent, when comparing eqs. (3.18) and (3.14), that 


for the given load W, and angular velocity w, dynamic 
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38 
bearing forces, P{ and P? are functions of: e€o,e,e,W and W. 


Now, for the given system parameters, if initial values 
for e, e€, Ww and w are known or assumed, the dynamic. bearing 
forces can be calculated using egs. (3.14) afterwhich the 
equations of motion (3.7) can be solved. Repeating the 
calculations in a time stepping manner, the solution for the 
shaft centre motion and dynamic oil-film forces is advanced 
in time. By non-dimenSionalizing eqs. (3.7) it can be shown 
[7] that, for constant angular velocity w, the shaft center 
motion iS governed by three independent non-dimensional 


parameters, namely: eo, e/C and g/Cw’. 


Considering the case when angular shaft velocity is 
time dependent, the equations of motion for the journal 
centre may be written 

P4+WcooswW+Mr-eQ*cos(Q-w)+Mr-eNMsin(Q-W) =MrC(e-W7e) 
(350099 
P$-Wsinw+Mr-eQ?sin(Q-wW)-Mr-eMcos(Q-W) =MrC(ewt+2e) 
These equations differ from eqs. (3.7) only by one 
additional term,- that iS an inertia force due to. shaft 


Circumferential acceleration. 


The Reynolds' equation of lubrication (3.8) is derived 
from the general Navier-Stokes equations [35] by considering 
several simplifying assumptions. One of the assumptions 
underlying Enuoe derivation is) thats “anertia ‘forces of 
lubricant are negligible. These inertia forces consist of 


fluid gravity, centrifugal forces acting in curved films and 
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acceleration of the fluid. This means that the lubrication 
equation ignores the inertia term due to unsteady velocity 
of the journal surface. During start-up operation, there is 
always tangential acceleration of shaft. The effect of this 
acceleration® tends to reduce the magnitude of oil-film 
pressure which, consequently, reduces instantaneous’ load 
Carrying capacity of journal bearing. This effect can be 
accounted for by including some of the terms originally 
dropped from the Navier-Stokes equations. This, however, 
would considerably increase difficulty of the problem and 
considerably augment the cost of the numerical solution. 
Therefore, in the present analysis oil-film forces are 
calculated assuming that the Reynolds' lubrication equation 


holds for transient journal angular velocity. 


3.4.3 Forcing function 


Equations (3.19) are solved using the fourth-order 
Runge-Kutta step procedure [37]. The dynamic film forces, Pj§ 
and P?, and their resultant in the vertical direction, are 
Calculated at each time step. A force, Fy, equal in 
magnitude and opposite in sign to this resultant is taken as 
the force (due to rotor imbalance) transmitted to foundation 
through the bearing. As with Model 1 (see Fig. 3.2 and eq. 


3.4), the forcing function for this model can be written 


Pet) =3P4 aeMue yi Cor 20)) 


a ae ae ae 


* Simplified analysis of this effect is presented in [35] 
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Pnsorder to vavoldametal tow metalasiraction at early 
Sragemeotm “Startoup operations “a-sroter is lifted up by 
pressurized oil supplied into journal bearing. This pressure 
is dropped when rotor angular velocity is high enough to 
secure sufficient load carrying capacity of bearing oil 
film. In view of the above, transient analysis of this model 


always starts with a rotor having some initial speed, w;. 


Thus, applying the reasoning discussed in connection with 
Modeiey sb 7(Sectuons eo 2a), rotor angular acceleration, 
velocity and acceleration are assumed to vary with time 


according to the formula: 
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4, COMPUTER PROGRAM DESCRIPTION 


In order to implement the reccurence formulae (2.23) in 
a computer program, components of generalized load vector 
(P;, Pz, and P3;) have to be evaluated at each time step, 


sccording to'eq. (2.16). 


Consider first Model 1 with the forcing function given 


by eq. (3.4). The following integration is required. 


ae 46 

Doar ocd t ko = Mr J. ¥m(t)at 
Th = z 

PR=f tF,(t)dt - Mri tym(t)dt (a1) 
T 


Deter (tt) at: a Mrf t?ym(t)dt 


The first integral in each of the above equations can be 
easily evaluated since F,(t) is a known function of time, as 
described in eqs. (3.3) and (3.5). The computer program 
employs a 4-point Gaussian quadrature integration scheme; 
results are stored in vector Q. In the process of deriving 
the finite time formulation (Chapter 2), shape functions N 
and their derivatives are determined, and vector q(t) is 
approximated by q(t)=[N]Jiq,} in eq. (2.9). The function 
Ym(t), aS a component of vector q(t), is assumed to vary in 
the time interval, +r, according to the same approximation. 


Simple integration yields 
PT=Q7-Mr (q?-qQ7) 
PP=Q3-Mr (rq?-qrt+tq?) (4.2) 
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The superscript "m" in eqs. (4.1) and (4.2) indicates that 
the above quantities represent only one component 
(corresponding to nodal variable ym) in each of the global 
(nx?) vectors. Substituting eq. (4.2) into the reccurence 


scheme (2.32) leads to the formulae: 


{q2} ois {Q,} {q,} 
fag? |e [Za ter) = aa Ga (4.3) 
{q2} Seed | ipa "TV {qy} 


which is actually employed in the computer program to 
advance transient solution in time. The subscript "m" 
indicates that the matrices Z,,; and 2Z;2 are "modified", by 
adding and subtracting from their appropriate components the 


terms given by eq. (4.2). 


The procedure just described is also used in the 
analysis of Model 2. In this case, however, two components 
(Fo, 8Fy) Of®the foreing? function given by ‘eq. (376) have “to 
be considered. Accordingly, two sets of equations’ must be 
integrated to evaluate components of generalized load vector 
P, As a result, vector Q has additional components and 


matrices Z,; and Z,2 are modified by additional terms. 


In the analysis of Model 3, matrices Z,; and 2;2 are 
modified in exactly the same way as in Model 1. However, the 


components of vector Q are obtained differently, since the 


* That is, eqs. (4.1) for the x-components and similar ones 
for y-components 
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termite, Gdintthenfiorcing lfiunctiom (eqn93e20): isunot given by 

an analytical expression. The procedure used in this case 

can be summarized as follows: 

in) Tueiueime gi ntervalienvecisrtdivided into  "j" ~‘sub- 
intervals 7. 

2. At each sub-interval 7, the dynamic force Fj is 
determined (according to the procedure outlined in 
Chapter 3). 

3. At the same time, two additional functions are 
generated, namely: 7F4 and 7?*Fj. 

4, Then, the "j" discrete values of these three 
functions: (ih. capers, "tks and vet*FQ) “are Hintegrated 
numerically over the interval +r to give components 
of vector Q. (The cubic spline method of numerical 


integration [38], is used in this process). 


4.1 Main Features of the Computer Program 


The basic features of the computer program are as 


follows: 


5. The program is written in FORTRAN language (double 
precision is used). 

6. It is general for any beam or plane frame models 
with horizontal and vertical beam members (can _ be 
easily extended for any structure). 

7. The program is general for any boundary conditions 


and constraints in the system. 
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Two different models of the Timoshenko beam finite 
element are used. 

The program can be used to determine system natural 
frequencies and eigenvectors. 

A static test can be performed, by computing static 
deflections (using global stiffness matrix) and 
comparing results with known exact solutions. 

The transient vibration problem of the models is 
solved using time stepping scheme based on finite 
time formulation. 

The program uses the following IMSL subroutines: 

a. EIGZS to determine eigenvalues/vectors 

b. LEQT2F as equation solver for static test 

Ga) LENVZE to invert matrix £245), 

The program employs 4-point Gaussian quadrature (for 
Model 1 & 2) and cubic spline (Model 3) as integrat- 
ing procedures to evaluate components of generalized 
load vector. 

The Runge-Kutta (4th order) method is used to solve 


journal centre equations of motion (Model 3). 
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4.2 The Computer Program Organization 


The simplified schematic diagram of the computer 


program structure is shown in Fig. 4.1 and Fig. 4.2. 


MAIN 


READ IN PARAMETERS 


(NBE;NCE;IE;IR;NL) 


CALCULATE MATRICES 


AND VECTORS DIMENSIONS 


if NCE>0 


| CALL reac | 


if NCE=0 


| CALL BEAM | 


NOTE! 1) Subroutines BEAM and FRAME are actually the main programs. 
The split MAIN-BEAM-FRAME is used to enable changing the 
dimensions of all global matrices and vectors, depending on a 


type and number of finite lements used. 


2) NCE - a number of elements in column. All other variables' 


Names are explained by comments cards in the program. 


3) The main segments and several subroutines of the computer 


program are listed in Appendix A-3. 


Fig. 4.1 Schematic diagram of the computer program (MAIN) 
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BEAM 
TAR 


READ IN PARAMETERS 
(NTS ;DT;ALF; BET;GAM; DEL) 


if IE=6 if IB=10 


READ IN [NG] ARRAY 


READ IN GLOBAL PARAMETERS 
(PRO;EMO; RHO; SHF;ECC;ETA) 


READ IN [NG] ARRA 


READ IN BEAM PARAMETERS 
(BML;ARE; SMA) 


if IE=6 if IE=10 


CALCULATE [EM];[EK] 


CALCULATE [EM]; [EK] 


ASSEMBLE [EM];[EK] INTO [GM]; [GK] 


MODIFY [GM] Pence 


INVERT F.T.E. SUBMATRIX [2,2] m 


INTEGRATE LOAD VECTOR 


EMPLOY RECCURENCE SCHEME 
RETURN 


Fig. 4.2 Schematic diagram of the computer program 
(subroutine BEAM) 
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5. SIMULATION RESULTS AND DISCUSSION 


In order. to make the discussion of the results 
(presented herein in graphical form) general, all the 
parameters and variables appearing on plots are 
dimensionless. Two parameters are especially important in 


the following discussion, namely: a and £, defined as: 


a=w,/Wo B=T1/To 


For the sake of simplicity, the dynamic response of each 
model is examined by analyzing displacements only at the 


Structure driving point (i.e. at the midspan of the beam). 


5.1 Parametric Studies - Model 1 


5.1.1 Preliminaries 


The dynamic analysis of a structural system is uSually 
preceded by determination of its natural frequencies. For 
Model 1, these frequencies depend on: (1) the supporting 
beam (i.e. its material, geometry and boundary conditions) 
and (2) the mass of the rotor. It is therefore of interest 
to know the effect of the rotor mass on the system natural 
frequencies. The natural frequencies of Model 1 are computed 
for different rotor to beam mass ratio with all other 
Parameters of the system fixed. The first three frequencies 
are plotted in Fig. 5.1. As expected, an increase of the 


mass ratio decreases the natural frequencies, except for 
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those of anti-symmetric modes, which are not affected. 
(Take, for instance, the second mode one indicated by "2" on 
the plot). Therefore, a fixed mass ratio 7=Mr/Mb=0.1, 
representative of large turbomachine support system, is 


chosen for further analysis. 


Fig. 5.2 shows build-up of the system vibration as a 
response to the exciting force generated by an accelerating 
unbalanced rotor starting from rest. Exactly the same 
response is plotted again, in Fig. 5.3, together with the 
"Dositive" envelope of maximum displacements. For clarity, 
only displacement envelopes are plotted on all further 
graphs. It should be noted that an envelope represents a 
series of points obtained from the numerical calculations. 
The shape of the curve, therefore, is not important since it 
depends on the approximation method chosen for connecting 


the associated data points. 


Fig. 5.4 shows the relationship between displacement 
envelope and time for three different lengths of the 
Supporting beam with, all other system parameters fixed. The 
plot is intended to show that the dynamic magnification 
factor (DLF) depends on the slenderness of the beam. In 
further analysis, only the beam _ shown in Fig. 3.3 (with 


Slenderness factor x=0.12) is considered. 
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5.1.2 Dynamic analysis 


Response of an ideal system without damping is examined 
first. An interesting case is shown in Fig. 5.5, which will 
be discussed in detail. This graph shows displacement 
envelopes versus time for three different values of the 
parameter 6B (i.e., for three different rotor acceleration 
times). The rotor operating speed is set to be greater than 


the first natural frequency of the system (a=1.2). 


Consider the case when $=45.0. Starting from zero, the 
displacement envelope increases slowly, but with a steadily 
increasing rate of change, following the increase of the 
rotor instantaneous speed, %. This rate of change reaches a 
maximum when the rotor speed approaches the fist critical 
frequency of the system (Q=w.), indicated by "e" on the 
curve. Passing this point, the displacement envelope 
increases further until gt reaches itS maximum at 
non-dimensional time t/7r,~29.2. (At this moment, the rotor 
instantaneous speed is still less than the operating speed 
w,, to be reached at time 45.0). After reaching its maximum 
value, the displacement envelope oscillates around a certain 
level with all the consecutive maxima slightly less than the 
Parst one. The spattern of ~this) oscillation ‘and its 
frequency, which depend on the rotor speed, stabilize after 
the rotor attains a constant operating angular velocity. For 
the fixed parameter a, as it is in the case being discussed, 


shortening the rotor acceleration time (f=25.0, 10.0) 
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increases itS acceleration rate. In comparing all three 
curves of this plot, it is evident that the maximum response 
amplitude and the level of the envelope oscillation are 
highly dependent upon the rotor acceleration rate. They 
become smaller when the acceleration rate increases, as does 
the difference between the maximum amplitude and the 
consecutive maxima of the envelope oscillation. It is also 
observed that the maximum amplitude does not occur at the 
cmidical mirequency;® andy» that a: shritein ats: position from 
that point also depends on the rotor acceleration rate. When 
the acceleration rate is high enough (as it happens to be 
for ~=10.0 on the graph), the maximum amplitude occurs some 
time after the rotor speed stabilizes at its operating 


level. 


Fig. 5.6 shows the results of the analysis carried out 
for varying the parameter a (a=0.5,0.8,1.2,1.5) and constant 
Baecep=30.00e dititkcsl cleampfirom this ‘graph thateithe: dynamic 
response and the maximum amplitude of vibration are 
dependent on the level of the rotor operating speed, and 
whether it is below or above the critical frequency of the 
System. For example, consider the cases for a=0.5 and 0.8. 
Both the level and the amplitude of the envelope oscillation 
increase when the rotor speed approaches the natural 
frequency of the system. The level of this oscillation could 
also be predicted through consideration of a steady-state 
response analysis. It should be noted that varying the 


parameter a, when 8 is constant, changes not only the level 
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of rotor operating angular velocity but also rate of its 
acceleration. It is, therefore, logical to compare the 
Sesuucsnm fom. thes Cases *ta=1.2) and® 195¢0fsFigim5.68to the 
previously discussed in detail results of Fig. 5.5. The 
results of Fig. 5.6 confirm once again that the maximum 
amplitude of vibration and the level of envelope oscillation 
decrease when the acceleration rate increases. Additionally, 
they show that the amplitude of the envelope oscillation 
becomes smaller for greater parameter a, i.e. when "the 
spread" between the rotor operating speed and the system 


Critical frequency increases. 


The analysis has proved so far that there is a 
pronounced effect of the rotor acceleration rate on the 
dynamic response of the system. To determine this effect 
more precisely, new dimensionless variables are introduced 
into the analysis, namely: the rotor acceleration rate 
ehrough atheloteritical frequency et, «and» thelshiftt of the 


maximum amplitude of vibration o, defined as: 
E=Q29/Wo’ F=Qm/Wo 


where: 9%, is the rotor acceleration rate at the moment’ the 
rotor instantaneous speed passes through the first natural 
frequency of the system (Q=w.), and 2, is the rotor 
instantaneous speed at the moment the response amplitude 


reaches itS maximum. 
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The effect of the rotor acceleration rate through the 
Critical frequency on the maximum amplitude of vibration is 
SHOWMelMeh ig. sha/earihe effect of ‘this sacceleration on ‘the 
shift of maximum amplitude from the critical frequency, is 
presented in Fig. 5.8. The first relationship (Fig. 5.7) 
demonstrates clearly that, for a fixed level of rotor 
Operating speed (i.e. for constant parameter a), the greater 
the acceleration rate through the critical frequency the 
smaller the maximum amplitude of the system response. It 
also shows that, for the same values of this acceleration 
rate, the maximum amplitude is reduced by increasing the 
deviation of the rotor operating speed from the critical 
frequency, i.e. for greater parameter a. Fig. 5.8 shows that 
the shift in location of the maximum amplitude from the 
Critical frequency increases for higher acceleration rates 
and/or for higher levels of the rotor operating speed. For 
the case when a=1.1, with the acceleration rate exceeding 
some critical value (&>=5.0°10°7?), the condition is met when 
the maximum amplitude occurs some time after the rotor speed 
has stabilized at its operating level. This is marked on the 
plot by a dotted line. (An example of response envelope, for 


Sichecase,) 2s shown in, Fig. 5.5,ecurvestor §6=10.0). 


The relationships. tpresentedy inheFigse.e 5icyeandynoe8 
include the effect of the rate of change of acceleration 
when the rotor instantaneous speed passes through the 
Critical frequency. To eliminate this effect, the analysis 


is repeated once again, assuming this time that the forcing 
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function is generated by an unbalanced rotor accelerating at 
a constant rate (%=constant). The results, for a fixed level 
of the rotor operating speed (a=1.2), are shown in Fig. 5.9. 
It is evident that an increase in the acceleration rate 
reduces the maximum amplitude of response and shifts it from 
the critical frequency. It is observed that, for very small 
values of the acceleration, the shift of the maximum 
amplitude (c=2,/wo) is slightly less than 1.0, which implies 
that the maximum amplitude occurs just before the 


instantaneous rotor speed reaches its critical frequency. 


The effect of damping, always present in structural 
systems, 1S now examined. Fig 5.10 shows the relationship 
between the displacement envelope and time, for a=0.8, 1.1 
and 1.2, $=30.0 and §=0.02, where § is a damping factor in 
the first mode of the system natural vibration. The effect 
of damping on system dynamic response can be clearly seen by 
comparing curves (with equal values of parameter a) shown in 
Fig. 5.10 and Fig. 5.6. Consider first the case a=0.8 (that 
is, when the rotor operating speed is less than the first 
Critical frequency). The response envelopes for the system 
with and without damping increase almost identically until 
they reach their maximum values. Then, for both cases, the 
envelope starts to oscillate about the same level which 
could be obtained from steady-state analysis. The effect of 
damping is to diminish the amplitude of this oscillation 
until it dies out completely. For the case with rotor 


Operating speed greater than critical frequency, the effect 
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of damping is much more pronounced. Comparing the results 
Bhovoeernetritcsers. 10Thend FigkSt6ésonce again (this time for 
a=1.2), it is noted that the maximum amplitude of vibration 
wow ereauced by 29%) Ubrom 15.5. to 1170), and the level of the 
envelope transient oscillation by 84% (from 12.5 to 2.0). 
For the same amount of damping in the system, the maximum 
amplitude and level of envelope transient oscillation 
depends on the rotor operating speed. This is clearly 


illusteatedtiby cunvesetfor éa=1elitandchre2)iginthrger5st10e 


Fig. 5.11 demonstrates the relationship between the 
maximum amplitude of vibration and the damping factor, for 
various levels of the rotor operating speed. The maximum 
amplitude of damped vibration is non-dimensionalized by 
dividing by the maximum amplitude of vibration of the 
undamped system subjected to the same excitation. It is 
observed that the maximum amplitude of vibration decreases 
with an increase in the amount of damping in the system. The 
effect of damping on reducing the response maximum amplitude 
becomes greater if the rotor operating speed is "closer" to 
the system natural frequency, that is for values of 


parameter a closer to 1.0. 


5.1.3 Concluding remarks 


The results obtained in this section are identical to 
hose neot etVuctor c& Elitvyintioh, andwthus: confirm thetabrlity 


of the model to predict known results accurately. Several 
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conclusions can be drawn from the results of numerical 


analysis of Model 1. These include: 


1. 


The maximum amplitude of vibration of low-tuned 
foundation structures Supporting rotating machinery 
occurs during the transients at start-up or 


shut-down operations. 


The response maximum amplitude is dependent on: (1) 
the rotor acceleration rate through the critical 
frequency of the system, and (2) the deviation of 
the rotor operating speed from the critical 
frequency. 

a. The greater this acceleration, the smaller’ the 

maximum amplitude. 
b. The greater the deviation of the rotor operating 


speed, the smaller the maximum amplitude. 


Mheremei:Si a yvshift. "1n ~ position = ot ~ the maximum 

amplitude with respect to the critical frequency, 

dependent on the same parameters. 

a. The greater the rotor acceleration rate through 
the critical frequency, the greater the shift. 

b. The greater the deviation of the rotor speed 
from the critical frequency, the greater the 


shift. 


The effect of structural damping is to subdue 


transient vibration to the level of steady-state 
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response, which depends on the rotor operating 
Speed. The maximum amplitude of vibration decreases 
with an increase of the damping in the system. The 
effect of damping on this reduction becomes greater 
if the deviation of the rotor operating speed from 


the critical frequency is smaller. 
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5.2 Parametric Studies - Model 2 


5.2.1 Preliminaries 


As in the case of Model 1, the dynamic analysis of this 
model is also preceded by the determination of system 
natural frequencies. For a simply supported beam, once its 
material and boundary conditions are defined, the natural 
frequencies spectrum depends on one non-dimensional 
parameter x. However, there is no such simple relationship 
for a framework. For a single bay portal frame (that is, for 
the simplest case of a framework), with clamped columns, the 
natural frequencies spectrum is a function of a group of 
several non-dimensional parameters (3:9; ]F. Therefore, 
parametric studies of Model 2 cannot be generalized as 
eaSily as they can be for simple beams. To simplify the 
discussion, a new non-dimensional parameter ("frame 
Bacanecer., vy), sidefined “aS (y=Ki/kK2-VlyAr, loa Le/b1, is 


brought into the analysis. 


One of the main distinctions between the vibration of a 
beam and a framework is that in the latter case, the 
longitudinal motion is coupled with flexural motion. 
Fig. 5.12 shows the distribution of the lowest four natural 
frequencies of the single bay portal frame (shown in 
Fig. 3.5), aS a function of the frame parameter y. The 
natural frequencies are obtained by considering the 


axial-flexural coupling in the model. The frequencies of 
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symmetric and asymmetric modes of vibration are indicated on 
the graph by "S" and "A", respectively. In order to simplify 
dynamic analysis of a framework, it iS a common practice to 
neglect the axial deformation of its members and to take 
only flexural motion into consideration. This approximation 
is applied to the frame model. The resulting natural 
frequencies of the system are plotted in Fig. 5.13. For easy 
comparison of results , the lowest three natural frequencies 
Ore rigs. S.t2 and.o.l3e8are: plotted together in Fig. 5.14. It 
is noted that as the frame parameter increases, the 
influence of axial deformation decreases, with flexural 
deflections eventually dominating the first three modes of 
vibration. However, a more detailed analysis indicates that 
the higher modes remain affected by axial vibration. As 
well, this approximation is not justified at all for small 
values of the frame parameter (clearly seen in Fig. 5.14). 
In view of the above, the axial-flexural coupling in the 


frame model is considered in the present analysis. 


Pig.@5, Seshowstthe effect Fot ether votoriito Ssupport 
frame mass ratio (yn=Mr/Mf£) on the determination of the 
lowest two natural frequencies of the system, for two 
different lengths of frame columns. This figure is intended 
to exemplify the reduction of natural frequencies’ for 
increasing ratio 7, and to demonstrate that this effect 
depends on the frame geometry and mode of vibration. The 
lowest four natural frequencies of the system chosen for the 


dynamic analysis (y=4.8) are plotted against the mass 
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rataoyi tam Fig.tios 16.0 Oncei againwva ‘pronounced; effect ‘of 
this mass ratio on the system natural frequencies can be 
seen. In-depth analysis of this effect would require 
Simultaneous examination of the system mode shapes of 
vibration, which is beyond the scope of this project. 
Further analysis was carried out for the fixed value of the 
mass! ratdo,s 7=0.:05;, representative OF turbomachinery 


foundation systems. 


5.2.2 Dynamic analysis 


The displacements at the location of the driving force 
(i.e. at the midspan of the beam) versus time are plotted in 
Fig. OralaLs Displacements in horizontal and vertical 
directions are marked on the graph by "X", and "Y", 
respectively. The displacements in both these directions are 
non-dimensionalized by dividing their absolute values by the 
Same quantity- y,. Where y, is the maximum static deflection 
at the midspan of the beam due to a vertical load equal to 
Mr-e-w,* applied at this point. The rotor operating speed is 
set to be beyond the fourth natural frequency of the system, 
i CeO=G)/Og=1.2.08nAS | bllustrated ein Fig. 5.12, the fourth 
natural frequency for this specific model (y=4.8) is of the 
Symmetric mode of vibration. The frame chosen’ for the 
dynamic analysis has a short stubby beam and relatively long 
and slender columps,e AS al@result of this’ geometric 
configuration and the specific frequency of the forcing 


function, the response of the system (at the driving point) 
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is dominated (in the y-direction) by flexural vibration. The 
motion of this point in the x-direction is more complex, due 
to superposition of the columns' flexural and the beam's 
axial vibration, which results in "beat" frequencies. The 
maximum displacement envelopes, for the same response, are 
Shown in Fig. 5.18. The instant when the instantaneous rotor 
Speed passes consecutive critical frequencies of the system 
is indicated on the curve by "e", There are four’ such 
poantspzacdescribedaby “lA "2S "ea "GAT ,cande "4S. wheres "Av 
and "S" stand for asymmetric and symmetric modes of 
vibration, respectively. It 1S observed that at the 
beginning of the transient period, the enevelope of maximum 
displacements in the vertical direction rises very slowly. 
For t/rTo<=2.0, the response is dominated by vibration in the 
x-direction. When the instantaneous rotor speed passes the 
third natural frequency, the vibration in the y-direction 
begamsatto adbualdeupremore wrapidly.» Eventually, ) @both: the 
maximum amplitude of vibration and level of envelope 
oscillation in the vertical direction become several times 
greater than in horizontal direction. Therefore it is 
logical to focus attention on the response envelope in the 


y-O1rection, forithis particular model. 


Fig. 5.19 shows the displacement envelopes versus time 
for a fixed level of rotor operating speed (a=1i.1) and 
Varioustemotor “acceleration times T,, (6=3.0,4.0,5.0).. The 
results, obtained for a fixed rotor acceleration time 


(B=5.0) and different levels of rotor operating speed 
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a= 1 cael. 2a aremlepresentedsrin® Fig. ©5720.¥eCareful 
examination of the curves shown in these two figures 
indicates that the patterns of the envelope oscillation are 
identical to those presented in Figs. 5.5 and 5.6 (for Model 
1). Evidently, the maximum amplitude of vibration and its 
shift from the critical frequency are both dependent on the 
level of rotor operating speed, a, and the rotor 
acceleration time, $8. It is observed that this dependency 
has exactly the same character as the one already discussed 
in detail in Section-5.1.2. It is obvious that the resulting 
general relationships between the maximum amplitude of 
vibration (and its shift) and the rotor acceleration rate 
througheathe | ycriticalefrequencyMof the systemiare; forythis 
model;Msimibarmito those shown in 3Figs. 5.7, 5.8 and 5.9 
(that is, for Model 1). Therefore, because they are "costly" 
(i.e., many runs of a computer program are required to 
generate data points), these general relationships are not 


plotted for Model 2. 


The effect of damping on the system response is shown 
in Fig. 5.21. The level of rotor operating speed and its 
acceleration times ace  faxedmeda-iem, «f=o.0)), ewitho se the 
displacement envelopes obtained for ideal (§=0.0) and damped 
(§=0.1) systems. The effect of damping is clearly seen and, 
Since it is identical to Model 1, no detailed discussion is 


presented. 
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An interesting case is shown in Fig. 5.22. Here, the 
rotor operating speed is set to pass the third natural 
Frequency of the system, which for this model is of 
asymmetric mode of vibration (see Fig. 5.12). As a result, 
the dynamic response of this specific model is dominated by 


VaDratiron in the horizontal ("xX") sdirection. 


5.2.3 Concluding remarks 


The following conclusions may be drawn from the results 


of the numerical analysis of Model 2: 


i. The proposed method for transient analysis proves to 
be equally suitable for frameworks as it is for 


Simple beam. 


2. The dynamic behavior of a framework, as a model of 
foundation structure, is far more complex than the 


response of a simple beam. 


3. Results of dynamic analysis of frameworks, even for 
Similar types of structures, cannot be generalized 
Since the dynamic response may differ qualitatively, 


depending on the specific model geometry. 


4, For a specific framework, as a model of a_ low-tuned 
Structure supporting rotating machinery, all the 
conclusions (i.e. regarding the maximum amplitude of 
vibration and the shift in its position, as well as 


the effect of damping on the system response) drawn 
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also valid. 
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Fig. 5.12 Effect of frame parameter on natural 
frequencies of a portal frame; axial 
deformation in beams considered. 
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frequencies of a portal frame; axial 
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5.3 Parametric Studies - Model 3 


5.3.1 Preliminaries 


At any given constant angular speed of an _ ideal 
(perfectly balanced) rotor, the journal centre is ina 
Stationary equilibrium and the total oil-film pressure force 
equals the static load on the bearing. However, when the 
rotor has an unbalance, the centre of the rotating journal 
is in motion describing a closed orbit about this 
equil@bniumayposition sc .Duenstonethis motion, additional 
pressures are set up in the lubricant film acting on the 
journal as dynamic forces over and above the static force. 
These dynamic forces, aS mentioned in Section 3.4.2, are a 
comphicatedhefunction Gof@megy, ene epory and V, where 
geen, LsRaG, Wee). tis describedndby eqs. 1.03015itand,/(3017)% 
This means that for specific rotor-bearing configuration and 
speed, the dynamic forces depend on both amplitude and 
velocity of the journal centre motion, i.e. on the journal 


centre orbit. 


FIG.) .5.2o Shows whirl -locisofjthe) journal ‘centre “for 
ep Ueemrand 9/Cw.=0525, being Styoical Sfor a rotor-bearing 
system of a relatively small high-speed turbogenerator'®. 
The orbits are obtained by integrating equations of motion 


(3.7) for two different values of the rotor unbalance 
'°Specific bearing system parameters used in the numerical 
calculations are taken from commercially operating turbo- 
generator (Wabamun Power Plant, Alberta, Canada), with the 
service speed 3600 RPM. 
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parameter, (e/C=0.1 and 0.2). The pronounced effect of rotor 
unbalance on the journal centre motion and, consequently, on 
o1l-film dynamic forces is obvious from this figure. In the 
dynamic analysis of previous models (that is with rigid 
bearings), the magnitude of the rotor unbalance was 
unimportant, Since it had no effect on the displacement 
envelopes! tri Thesuresultsheshowns tin Figers=5e23inindicate, 
however, that this magnitude will play a very important role 


in the analysis of the present model. 


Under dynamic conditions the bearing oil film behaves 
aS a spring-damping system. By employing a conventional 
linear analysis'? of the journal bearing, it can be_- shown 
that the oil-film "stiffness" and "damping" properties 
depend on the bearing geometry, the lubricant viscosity, the 
Static load on the bearing and, most importantly, on the 
rotor speed. As a result of these dynamic properties, the 
oil-film plays a dominant role in attenuating or amplifying 
the excitation force due to the rotor unbalance. These 
effects are very complex in nature, involving a great number 
of independent parameters, and their detailed discussion is 


beyond the scope of this project'®. 


mee ei ei 


‘tNote that the displacement envelopes are normalized with 
respect to static deflections due to centrifugal force 
Mr-ew;. 

'2Linear analysis of the journal-bearing system is very well 
documented. See, for example, [24-28]. 

‘23An interested reader is referred to specialized papers 
concerned with stability of the rotor-bearing system. For 
example [33-34]. 
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5.3.2 Dynamic analysis 

Fig. 5.24 shows the relationship between displacement 
envelope and time for the model with: (1) journal bearing, 
and (2) rigid bearing. Except for parameters associated with 
the bearing, all other system parameters are identical for 
both cases. The curves illustrate the results of the 
analysis carried out for the specific rotor-bearing system 
with the operating speed of 3600 RPM. The first natural 
frequency of the supporting beam selected for this analysis 
Gives a=w,/wo=1.2. The rotor acceleration time is set so 
thate the? ,.5sotor service speed is reached at the 
non-dimensional time equal to 20.0, (f$=20.0). Due to the 
reasons explained in Section 3.4.3, transient analysis. of 
these models starts with a rotor having an initial speed wo, 
(set to be 1000 RPM in this case). Therefore, transient 
analysis has to _ be preceded by initial calculations (time 
t=0) to establish the journal centre position and velocity 
for model (1), and dynamic forces transmitted to the 
Supporting beam for both models. Applying these initial 
forces to the beam, being at time t=0 in static equilibrium, 
results in shock excitation which causes some irregularities 
in the system dynamic response during a_ short initial 
period. Therefore, displacement envelopes shown in Fig. 5.24 


do not start at t=0. 


Comparison of the curves indicates that the dynamic 
response of the system has practically the same general 


character for both models considered. However, magnitudes of 
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the maximum amplitude of vibration and levels of envelope 
oscillation differ considerably. It should be pointed out 
that the absolute value of the rotor unbalance is identical 
in both models, and consequently, the response is normalized 
with respect to the same value. Hence, it is concluded that 
the force due to the rotor unbalance is magnified by the 
bearing oil film. It is also noted that the magnitude of the 
envelope oscillation after the rotor speed has stabilized at 
the operating level is almost the same for both cases. This 
suggests that the steady-state bearing transmissibility is 
close to 1.0 and, consequently, that the force magnification 


occurs during the transient period. 


Fig. 5.25 shows transient orbits of the journal centre, 
for this specific case, following the rotor coming up from 
the initial to the operating speed. High rotor acceleration 
rate at the beginning of the motion causes rapid increase in 
the rotor speed and, accordingly, fast growth of the orbit 
amplitude. With the steadily decreasing acceleration rate 
these changes become smaller for each consecutive cycle of 
journal motion and, eventually, the journal centre describes 
the steady-state orbit. From this illustration the effect of 
changing rotor speed on journal centre motion is apparent. 
From the remarks made earlier it follows that the dynamic 
bearing force and oil-film flexibility and damping have to 


change accordingly. 
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For the particular case studied the change of bearing 
dynamic forces with time is demonstrated in Fig. 5.26. Both 
forces (i.e. journal bearing and rigid bearing forces) are 
normalized with respect to the same value of maximum 
centrifugal force due to the rotor unbalance (Mr-ew?). Both 
nonlinearity of the oil-film force, and the change of 
transmissibility with time is noted. To better visualize the 
dynamic transmissibility of the Onli } a force 
magnification factor is defined as a peak-to-peak magnitude 
ratio of journal bearing force to rigid bearing force. This 
magnification factor versus time is plotted in Fig. 5.27. 
The relationship manifests ED agreement with the 
conclusions drawn from the analysis of the system response 


presented in Fig. 5.24. 


The effect of damping on the response is demonstrated 
in Fig. 5.28. The displacement envelopes for the model with 
the journal bearing and with the rigid bearing are _ plotted 
both for ideal and damped systems. It is noted that the 
levels of the damped transient envelope oscillation for both 
these models lie close to each other. This can be expected, 
Since the steady-state oil-film force magnification factor 


is less than 1.1, as shown in Fig. 5.27. 


To appreciate the effect of the rotor unbalance on the 
system response, the analysis was repeated for the rotor 
unbalance parameter e/C=0.1 with all other system parameters 


unchanged. The transient motion of the journal centre is 
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shown in Fig. 5.29. Comparing the results presented in Figs. 
5.29 and 5.25, it becomes obvious that transient orbits of 
the journal centre are highly dependent on the magnitude of 
the rotor unbalance. The smaller this unbalance, the smaller 
the amplitude of the journal orbit and, obviously, the 
smaller the magnitude of the dynamic bearing forces. 
However, no conclusions regarding the dependency of the 
bearing transmissibility on the rotor unbalance can be drawn 
from this figure. To do this, the instantaneous peak-to-peak 
force magnification factor is calculated. The results are 
Snownim Dief 1 okat 5-330 .codtloasenelearsesthatn for, Sthejagiven 
rotor-bearing system parameters, the Oi1lnfidm 
transmissibility depends on not only the rotor instantaneous 
speed but also on the rotor unbalance. Moreover, there seems 
to be no unique relationship between the transmissibility 
and the magnitude of the oil-film dynamic force. Finally, 
the system response envelope versus time, for e/C=0.1 and 
Unepeare plotted in Fig. 5.31. "ltrs observed that, the 
maximum amplitude of vibration and the level of envelope 
oscillation are both higher for the smaller rotor unbalance. 
This does not come as a Surprise, remembering that the 
response is normalized'* with respect to two different 
values, and since the transmissibility is higher for e/C=0.1 
as illustrated in Fig. 5.30. It should be stressed however 
that the absolute system vibration levels decrease with a 
FeOuiece VON I Natniearoror unbalance, 


'4The curves represent the dynamic load factor as shown in 
Fingea S125 
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Fig. 5.32 shows the results of the analysis carried out 
for varying the rotor acceleration time parameter B (8=20.0, 
30.0, 40.0) and all other system parameters fixed. The 
results illustrate that the maximum amplitude of vibration 
and its shift from the system critical frequency are both 
dependent on the rotor acceleration time, that is also on 
the rotor acceleration rate through the critical frequency. 
It is noted that the character of this relationship is 
identical to that previously discussed in connection with 
the lanalysis tofeModel?lwande2andSee sFigs poss s5onando'5.19). 
This implies that the conclusions of Section 5.1.3 regarding 
the dependency of the system response on the rotor 
acceleration rate through the critical frequency and on the 
rotor operating speed are, in general, irrespective of the 


type of bearing in the system. 


The analysis carried out so far dealt with the 
bearing-rotor system parameters typical tesor ouaek ‘smat 
turbogenerator. Larger unitS operate generally at lower 
speeds and with higher static loads on bearings. Both these 
parameters, as discussed in Section 3.4.2, have an essential 
effect on the dynamic performance of the journal bearing. It 
is therefore of interest to investigate the dynamic response 
of the system with considerably different parameters. 
Consequently, further discussion is concerned with some of 
the results obtained for eo=0.7 and g/Cw?=1.0, being typical 
of* a medium-sized turbogenerator rotor [7]. The rotor 


operating speed for this model is 1800 RPM. 
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Fig. 5.33 shows a family of steady-state journal centre 
orbits obtained for different values of rotor unbalance. It 
is noted, comparing the curves shown in Figs. 5.33 and 5.23 
(for the corresponding values of rotor unbalance), that 
non-linear effects are more pronounced for greater €o. 
Transient orbits of the journal centre for the _ rotor 
accelerating from the initial speed of 500 RPM and with the 
rotor unbalance parameter e/C=0.4 is presented in Fig. 5.34. 
The effect of changing rotor speed on the orbit amplitude is 
well illustrated. The change in position of the 
"instantaneous" steady-state journal centre can be easily 
traced as it "moves" upwards following an increase in the 
rotor speed and, consequently, in the bearing load carrying 


Gapacrtys 


The dynamic response of the system to_ unbalance 
excitation, for the models with journal bearings and with 
Ploideebearings, as) shown. an sFigs 5235. The rotor speed 
parameter and acceleration time parameter are a=1.2 and 
B=20.0, respectively. While the general character of the 
System response remains the same for both models, the 
maximum amplitude of vibration, the level of envelope 
oscillation, and the amplitude of this oscillation are 
greatly increased for the model with journal bearings. This 
implies that the journal bearing transmissibility, for this 
specific case, is much higher than 1.0. More importantly, 
the increased magnitude of the envelope oscillation infers 


that the transmissibility remains high even after the rotor 
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94 
speed has stabilized at its operating level. 


Fg. 30 illustrates the relationship between 
displacement envelope and time for the same rotor-bearing 
system and a different supporting beam (a=0.8). The results 
presented are obtained for the rotor unbalance parameter 
e/c=0.1 and 0.4 and for the models with journal bearings and 
with rigid bearings. The curves indicate that the system 
response is typical for high-tuned supporting structures 
Subjected to unbalanced excitation, as discussed in detail 
in Section 5.1.2. They also imply that the transmissibility 


is higher for the smaller rotor unbalance. 


Fig. 5.37 shows the peak-to-peak force magnification 
factor plotted against time, for e/c=0.1 and 0.4. These 
relationships fully confirm the predictions deduced from the 


analysis of the results presented in Figs. 5.35 and 5.36. 


5.3.3 Concluding remarks 


The following major conclusions are drawn from the 


results of the numerical analysis of Model 3: 


1. The proposed method of solution for the transient 
response analysis of structures supporting rotating 
machinery is readily adaptable for complex models of 


Phe forcing ulupet ion. 


2° w Lhe dynamic response of the system is highly 


dependent on the rotor-bearing system parameters. 
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extremely OEE Loan’ 


independent parameters 


Dynamic properties of journal bearing depend on 


result 


overall trends in the system 
of the system parameters 
example, lubricant viscosity, 


clearance et cetera, is 


due to a large number of 


involved. 


the 


region of bearing operation on design maps (€o, 
g/Cw2) and they change during transient period with 
rotor speed. 

The magnitude of the rotor unbalance has a 
Significant effect on the system response. This 
effect is dependent on a specific rotor-bearing 
System and varies with the instantaneous’ rotor 
speed. 

For any given system parameters, the dynamic 
response is dependent on the level of the _ rotor 


operating speed and 


the system critical 


regarding 
damping on the 


analysis of Model 1 


system 


its acceleration rate through 


frequency. The conclusions 


this dependency, as well as the effect of 


response, drawn from _ the 


(See Section 5.1.3) are general 


regardless of a type of bearing in the system. 
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Fig.25.9% Transient whirl orbit of journal centre 
due to unsteady rotor angular velocity. 
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6. CONCLUSIONS AND RECOMMENDATIONS 


A numerical investigation was carried out mainly to 
verify the proposed method of solution and to establish its 
usefulness for the specific type of problem. Both the 
results and the simplicity in obtaining them demonstrate 
that the method of solution based on finite time formulation 
is very convenient for the transient analysis of machine- 
foundation interacting systems. The study shows that’ the 
method is suitable for any system regardless of structure 


enay,or forcing function complexity. 


General conclusions, drawn from the numerical analysis 
of each of the models considered, regarding dynamic behavior 
of the system are given in the sections immediately follow- 
ing the discussion of the results. Perhaps the most import- 
ant conclusion being that the maximum amplitude of vibration 
of low-tuned structures supporting rotating machinery occurs 
during the transient period of the rotor coming up to speed 
and that it is highly dependent on the rotor acceleration 
rate through the critical frequency of the system. This 
maximum amplitude cannot be predicted by the classical 
Steady-state analysis. Therefore, the transient response 
analysis, as the most inclusive approach to the system dyna- 
mic analysis, should be employed in the present-day design 


practices of low-tuned foundations. 


Models chosen for this study were simple ones. They 


proved to be adequate for a qualitative analysis concerned 
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mainly with the determination of general relationships and 
trends in the dynamic behavior of the system. It should be 
stressed that more accurate models would not alter the gene- 
ral conclusions drawn from this study. However, for a quan- 
titative analysis of the system response, the use of much 
more elaborate models would be required. That means real- 
istic modelling of three dimensional foundations, multiple 
bearings flexible rotor systems, bearing pedestals, seals, 
et cetera. While such a comprehensive approach to the 
problem is very desirable and possible with the analytical 
and technological tools available today, it is not always 
practical. The analysis would be extremely difficult due to 
the further increase in the number of independent parameters 
involved and because of a drastic increase in the size of 
the problem followed immediately by many numerical 
difficulties. Moreover, in many instances the cost of 
transient analysis, involving the repetition of tedious 
calculations for a great number of time stepS may simply 
prohibit the use of rigorous analysis of the system. As a 
result there will always be a need for intermediate models 
and methods allowing for a simplified analysis of the 


system. 


It is Suggested that the method of solution presented 
in this study constituteS a very encouraging basis’ for 
further development. It would be interesting, for example, 
to determine effects of changing bearing geometry and 


lubrication system parameters on the system transient 
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response. The method can also be used to obtain information 
on system stability, response to impact loading, and so on. 
Before such efforts are undertaken, thorough studies are 
required to determine the numerical accuracy and stability 


of the method. 
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APPENDIX A-1: Nomenclature 


denotes column vector (nx7) 

denotes square matrix (nxn) 

denotes transpose of square matrix (nxn) 
denotes reciprocal of square matrix (nxn) 
cross-sectional area 

constants 

Subscripts, refer to beam 

damping matrix 

radial clearence in journal-bearing system 
Subscripts, refer to column 

modulus of elasticity 

rotor mass eccentricity 

eccentricity of journal centre 

subscript, refers to frame 

external load vector 

Shear modulus 

gravitational constant 

oil-film thickness, h=C(1+ecos@) 

identity matrix 

second moment of cross-sectional area 
CounternGr=152)3;,.00) 

stifness matrix 

Shear coefficient 

cross-sectional radius of gyration, k=/I/A 


length (beam/column/bearing) 
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mass matrix 

rotor mass 

Shape (interpolation) functions 

oil-film pressure 

generalized forces 

centrifugal force due to rotor imbalance 
Oil-film dynamic forces 

nodal displacement, velocity and acceleration 
modified Sommerfeld Number 

time 

time interval (time step) 

rotor acceleration time 

nodal displacement (axial) in beam element 
Shear force 

nodal displacement (transversal) in beam element 
rotor weight 

displacement (horizontal) at beam midspan 
displacement (vertical) at beam midspan 


matrix of finite time formulation 


rotor speed parameter (frequency ratio), a=w,/wo 
rotor acceleration time parameter, B=T;/To 
frame geometry parameter, y=/1,A2/12A,'L2/L, 


journal eccentricity ratio, e=e;/C 
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Eo eccentricity of static journal centre 

§ damping factor 

n rotor to support mass ratio 

6 cross-sectional rotation (in Thimoshenko beam) 

6 circumferential coordinate (in journal-bearing) 

0; ,0> limits of oil-film positive pressure 

K beam slenderness factor, k=k/L 

mM viscosity 

v Poisson's ratio 

3 rotor acceleration (through the critical frequency) 


Parameter, £=2 /wo? 


T Constant. t=3. 1415927... 

p density 

0 shift of maximum amplitude of vibration, o=Q2,,/wo 
T time interval (time step) 

7 time sub-interval 

16 system first natural period 

p bending slope 

x shear coefficient (on graphs only) 

YW shear slope (in Timoshenko beam element) 

vy) attitude angle (in journal-bearing system) 
W,W, rotor operating (service) speed 

Wr rotor operating speed (on graphs only) 

Wj rotor initial speed 

Wo System first natural frequency 


Wn Natural frequency of the n-th mode of vibration 
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rotor angular travel 

instantaneous rotor speed and acceleration rate 
rotor acceleration rate through the critical freq. 
instantaneous rotor speed at the instant the 


response amplitude reaches its maximum 


APPENDIX A-2: Timoshenko beam elements' matrices 


Simple Timoshenko Beam Element 
(Davas, atvaly-[1s)) 


(a) element stiffness "matrix: 


(EK |] =a 
k, 0 0 
Symmetric ke —ks; 
Ky 
where: 
@ = 12EI /GKAL? 


@. =3bP/(1+¢)L> 


k, = AE/La 
k 2 = 12 
k 3 = 6L 


Ky = (4+@)L? 


ks = (2-¢)L? 
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(b) element mass matrix 


M3 My 0 Ms —-M65 
M7 0 Me¢ Mg 
[EM] = £8 
m, 0 0 
symmetric m3 My 
M7 
where: 


B = pAL/(1+¢)? 


We zal Al 


ples fei ssye 

Wey Elie 96% 

Mz-= 13/35+7¢/10+¢7/3t67/5 

m, = L(11/210+116/120+¢7/24+7/10-y¢/2) 

ms = 9/70+3¢/10+o7/6-6y/5 

me = L(13/420+3¢/40+67/24-y/10+7¢/2) 

m, = L*(1/105+¢/60+$2/120+2y7/15+7¢/6+7¢"/3) 
me = L2(-1/140-¢/60-¢7/120-7/30-y¢/6+7¢? /6) 
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Complex Timoshenko Beam Element 


(Akella & Craggs, 


(a) element stiffness matrix: 


LEK] = 


where: 


0 0 0 0 
Ke kisaees kK 5 Ks 
Kem ke Kg 
Kapaa eae 
Ky 
symmetric 
GKA/L 
E/GK 
6/5 
375 
L/10A 
Teal OT 
3L*/10+6EI/5GKA 
L?/20A 


Ky 


aie) 


0 0 0 
k3 ky ks 
kp ky kis 
“kolo yee kas 
Kio ekki 
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=ky Uk? Wks 
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kp = L*/20I1+EL/10GKA 

ky = 3L7/10-6E1/5GKA 

Kao t= 720 lee / 1 0GKA 

kai = 4L°/30A- 

Ki2 = be/120TA 
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(b) element mass matrix: 


ny= 0 0 0 0 M2 0 0 0 0 


Mm; Neves Mg 0 Mz; —-Mg Mg Mes 
Mio —M;; M12 O Meg ~—M,3 Miy Mi5 
Mis —M1,7 O —Mg Miyg ~M,g ~M17 
Mig 0 Mp ~“Mi5 M17 M20 
[EM] = o 
mM; 0 0 0 0 
m3; —Myg Ms M6> 
symmetric Ne peah ae ep re oe Nchae: 
Mi5 Ki7 
M49 
where: 


0) = pAL 
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1/3 

1/6 

1365 

77280 

11L/210A 

L?/2401 

9/70 

11L/280 

13L/420A 
17L7/1260+131/35A 
57L7/5040A 
11L7/100801+11L/210A 
L?/90-91/70A 

Le a2A 
11L7/100801-13L/420A 
E27 4.05R" 

9L°/10080AI 

L?/140A? 

i 0080 N44 77/7 105A 
Ey 100601 *-L7/ 140A1 
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APPENDIX A-3: Computer program listing 


The following is the computer program listing (in 
FORTRAN code) consisting of the main segments (TRANSFB, BEAM 
and FRAME) and several subroutines. The schematic diagrams 
of TRANSFB and BEAM are shown in Figs. 4.1 and 4.2. There 
are considerable number of comments cards included in the 
program. It is hoped that these comments cards together with 
the above mentioned diagrams and the theoretical discussion 
presented in Chapters 2, 3 and 4 should make the program 


self-explanatory. 
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0002 
0003 
0004 
0005 
0006 


0007 


0008 
0009 
0010 
0011 
0012 
0013 
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program 


name: TRANSFB 


FEFEEEFLE EE EFTE EEE EE EELEEE LE LEEFLFEEFEFEPELE PETE ETH 


++ttteepeetetetetete te ttt 


OF SIMPLE BEAM OR PORTAL FRAME AS MODELS OF 
LOW-TUNED STRUCTURE SUPPORTING AN UNBALANCED 


A RECCURENCE SCHEME BASED ON FINITE TIME 


(SIMPLE OR COMPLEX) IS USED TO APPROXIMATE 
ELASTIC AND INERTIAL PROPERTIES OF STRUCTURE 


NATURAL FREQUENCIES AND/OR EIGENVECTOR 


TRANSIENT RESPONSE 


ACCELERATING ROTOR 


ale ke he 
PAS a PAS 


DIRECT INTEGRATION 
OF SYSTEM EQUATION OF MOTION USING 


ELEMENT FORMULATION 


1. 0s ' 
CLS ache RAS 


FINITE TIMOSHENKO BEAM ELEMENT 


1. (ns ‘, 
KR KW 


DETERMINATION OF SYSTEM 


t++ettetert ete tte t ttt 


AND STATIC TEST 


FEFFLFFEEFEFEEFEFEFEFELPEHEFPEEFEFEAPEPEE PEPE +++ TH+ H+ t+ 


REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 


INTEGER 


NE 
NBE 
NCE 
IE 
IR 
NL 
NG 


M(50,50) ,K(50,50) ,C (50,50) ,DD(50) 
AA (1275) ,BB(1275) 

E(10, 10) 

Z(150,150) ,V (150, 150) 

WV0 (150) ,WV1 (150) ,WV2(150) ,Q(150) 
WK (23000) 


NG (20, 10) 


total number of elements in the structure 
number of elements in the beam 

number of elements in each column 

number of D.O.FsS per element 

number of D.O.Fs of constrained structure 
nodal no. of applied load; vertical comp. 
array; elemments' nodal numbers in global 
numbering system; considering constraints 
and boundary conditions of the system 


READ (5, 1000) NBE,NCE,IE,IR,NL 


BS SR 


NE=NBE+NCE*2 
ISM=IR* (IR+1) /2.0+0.5 
MWK=IT**2+3%*I1T 


IF( NCE 
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53 
54 
35 
56 
57 
58 
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69 
70 
rie 
72 
as 
74 
ie 
76 
ie 
78 
a9 
80 
81 
82 
83 
84 
85 
86 
87 
88 
89 
90 
91 
SZ 
93 
94 
95 
96 
oF 
98 
og 
100 
101 
102 
103 
104 


0014 


0015 


0016 
0017 
0018 


0001 


0002 
0003 
0004 
0005 
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0007 
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0011 
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0014 


0015 
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CALL BEAM(NBE,NE,IE,IR,IT,NL,ISM,MWK,M,K,C,E,AA, 


+ BB,DD,WK,NG,Z,V,Q,WVO,WV1,WV2) 
10 CALL FRAME (NBE,NCE,NE,IE,IR,IT,NL,ISM,MWK,M,K,C, 
+ E,AA,BB,DD,WK,NG,Z,V,Q,WVO,WV1,WV2) 
1000 FORMAT (516) 
STOP 
END 
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The following subroutine is actually the main pro- 
gram. The split MAIN-BEAM was introduced to enable 
Changing dimensions of global matrices and vectors 
depending on a number of D.O.F. in the system. The 
BEAM subroutine differs slightly for each model 
considered. The version for Model 3 is as follows: 


SUBROUTINE BEAM(NBE,NE,IE,IR,IT,NL,ISM,MWK,M,K,C, 
+ E,AA,BB,DD,WK,NG,Z,V,Q,WVO,WVi,WV2) 


REAL*8 M(CIR,IR) ,K(IR,IR) ,C(IR,IR) ,ECIE, IE) 
REAE*Ss Z(t, PT) eVGT, IM)sBB Go) .,ce(so) 

REAL*8 WVOCIT) ,WV1 (IT) ,WV2(IT) ,QCIT) , WK (MWK) 
REAL*8 ALP,BET,GAM,TO,DT,TA,T,PI,OMO 

REAL*8 BML,EL,EMO,SMO,ARE,SMA,RHO,PRO,GKA,SHF 
REAL*8 MR,ECC,OMI,OMR,OMT,OMS,OMA,F1,F2,F3 


REAL*8 UC,WT,MC,ME,PAR,EPS,PSI,VEP,VPS,DX,PYI 
REAL*8 TX(30) ,PX(30) ,PY(30) ,VIS,BL,BR,RC 
REAL LAST 


INTEGER NG(NE, IE) 


— temporary element (mass or stiffness) matrix 
—- consistent mass matrix of constrained system 
consistent stiffness matrix of const. system 
- system damping matrix (C = ALF*M + BET*K) 


QA Ss 
| 


READ (5,1001)NDT,DT,ALP, BET, GAM 


NDT —- number of time steps 

DIB== SCIinersEeep 

ALP- constant in Rayleigh's proportional damping 
BET - constant in Rayleigh's proportional damping 
GAM - rotor acceleration time parameter 


WRITE (6, 1107) 
WRITE (6, 1114) 


If simple Timoshenko beam (6 D.O.F.) used goto 453 


IF(IE .EQ. 6) GO TO 453 
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106 
107 
108 
109 
110 
114 
pe 
Piss 
114 
TLS 
116 
Rt 
118 
wtS 
120 
WA 
TZ 
123 
124 
T25 
126 
127 
128 
29 
TG 
ta 
132 
$33 
134 
135 
136 
Foy 
138 
To 
140 
141 
142 
143 
144 
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148 
149 
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Bort 
(52Z 
13 
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0016 
0017 
0018 
0019 
0020 
0021 
0022 


0023 
0024 
0025 
0026 
0027 
0028 
0029 
0030 


0031 
0032 
0033 
0034 


0035 
0036 
0037 
0038 
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0040 
0041 
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453 
454 


55:1 


552 


tz 


WRITE (6, 1115) 

GO TO 454 

WRITE (6, 1116) 

CONTINUE 

WRITE (6, 1100) NE, NBE, NCE 
WRITE (6.41111) IE,IR,IT,NL 
WRITE (6, 1106) 


Read in NG array, storing element nodal numbers in 
global numbering system. If complex Timoshenko 
beam finite element (10 D.O.F.) used go to 551. 


TFCIE .EO. 10) "GO TO 551 

READ (5, 1000) ((NG(1,J) ,J=1,IE) ,I=1,NE) 

WRITE (6,1101) ((NG(1,J) ,J=1,IE) , 1=1,NE) 
GO TO 552 

READ (5,2000) ((NG(I,J) ,J=1, IE) ,I=1,NE) 

WRITE (6,2101) ((NG(1,J) ,J=1, IE) ,1=1,NE) 
CONTINUE 

WRITE (6, 1105) 


Read in beam material constants and shear factor 
PRO - Poisson's ratio 

EMO —- elasticity modulus 

RHO - density 

SHF - cross-sectional shape (shear) cofficient 


READ (5, 1002) PRO, EMO, RHO, SHF 
PI=3.141592654D+00 
SMO=0.5*EMO/ (1.0+PRO) 
RHO=RHO/ 386. 16 


Read in beam geometry parameters & mass of a rotor 
BML —- beam total length 

ARE - cross-sectional area 

SMA - second moment of cross-section area 

MR - rotor mass 


READ (5, 1002) BML, ARE,SMA,MR 

EL=BML/NBE 

WRITE (6, 1102) PRO, EMO, RHO,SHF,BML,ARE,SMA 
GKA=SMO*xSHF*ARE 


Initialize global M and K matrices. 


DO 10 I=1,IR 
DO 10 J=1,IR 
M(I,J)=0.0 

K (1,0) =020 


Call routine to evaluate complex (TM544) element 
stiffness matrix. If simple element used goto 661 
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IF(IE .EQ. 6) GO TO 661 

CALL KCTMBM (EL,GKA,ARE,EMO,SMA,E) 
GO TO 662 

CALL KSTMBM (EL,GKA,ARE,EMO,SMA,E) 
CONTINUE 


Call assemblage routine (looping for total number 
of elements in a beam) to get system stiff. matrix 


DO 12 NK=1,NBE 
CALL ASSZW(E,NK,IE,NE,IR,K,NG) 


Call routine to evaluate complex (TM544) element 
mass matrix. If Simpie element used go to 663. 


IF(IE .EQ. 6) GO TO 663 

CALL MCTMBM (EL, ARE,SMA,RHO,E) 

GO TO 664 

CALL MSTMBM (EL, GKA,ARE,EMO,SMA,RHO,E) 
CONTINUE 


Assemble elements' mass matrices into system's one 


DO 14 NK=1,NBE 
CALL ASSZW(E,NK, IE,NE,IR,M,NG) 


Read in journal-bearing system parameters 
Vis > — 2UubrPacant viscosicy 


BL - length of a bearing 

BR - radius of a bearing 

RC - radial clearence 

UC - normalized rotor unbalance (e/C) 


READ (5, 1002) VIS,BL,BR,RC,UC 
ECC=RC*UC 

RBM=MR/ (BML*ARE*RHO) 
WRITE(6,1122)MR,ECC,RBM 


Read in parameters 

OMO - system first natural frequency 
OMR - rotor oparating (service) speed 
OMIP=? rotor, anitiall, (speed 

HOM - system highest natural frequency 


READ (5, 1002)OMO0,OMR,OMI,HOM 
OMR=OMR*PI/30.0D+00 

DEL=OMR/OMO 

YO= (MR*ECC*OMR**2*BML**3) / (48 .0*xEMO*xSMA) 
TO= (PI*2.0D+00) /OMO 

TN=2.0%*PI/HOM 

TA=GAM*TO 
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ZET=DT/TN 


Read in number of time steps (IP) for calculation 
of oil-film dynamic forces and initial conditions 
IP - number of time sub-intervals 

EPS —- journal centre eccentricity 

PSI - attitude angle 

VEP — journal centre radial velocity 

VPS - journal centre tangential velocity 

PY SueVerercal Component Of initial dynamic force 


READ (5,1001) IP,EPS,PSI,VEP,VPS,PYI 
WT=MR*386.16D+00 

MC=MR*RC 

ME=MC*UC 

VIS=V1IS/6.894757D+03 
PAR=VIS*BL**3*BR/ (RC**2*2.OD+00) 
RSP=OMR*60.0/ (2.0*PI) 

SMN=PAR*OMR/ (WT*«2.0) 

STP=OMR* (RC/386.16) **0.5 
NDT2=TA/DT+1.5 

DT=TA/ (NDT2*1.0D+00) 
NDT=NDT2+TA/DT+1.5 

DX=DT/ (IP*1.0D+00) 

IP=IP+1 

WRITE (6, 1108) DT, NDT, ALP, BET, GAM, DEL 
WRITE (6, 1123)OMO0,TO,HOM,TN,ZET,OMR,TA 
WRITE (6, 1125) OMI 
OMI=OMI*PI/30.0D+00 


Evaluate system proportional damping matrix from 
global mass and stiffness matrices 


DO 30 I=1,IR 
DO 30 J=1,IR 
C(1I,J)=ALP*M (I,J) +BET*K (I,J) 
CONTINUE 


Call subroutines to create submatrices of a finite 
time formulation, and then modify them to include 
effect of rotor acceleration in its linear motion 
(due to beam vibration) on forcing function. 


CANE EI ZIZCIReT? DTM. CeKe2) 
Z(NL,NL+IR)=Z(NL,NL+IR)+MR 
Z(NLt+IR, NL) =Z(NL+IR,NL)-MR 

Z(NL+IR, NL+IR) =Z(NL+IR,NL+IR)+MR*DT 
Z(NL+2*IR, NL) =Z(NL+2*IR,NL) -MR*DT 

Z(NL+2*IR, NLtIR) =Z (NL+2*IR, NL+IR) +0.8*MR*DT**2 
Z(NL+2*IR, NL+2*IR) =Z (NL+2*IR, NL+2*IR) +MR*DT**3/60. 
IDGT=0 
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Invert modified finite time formulation matrix Z12 


CALE UENVERI( ZIT, TPA eT DCT, WK, TER) 
DO:25) 21 =r 

DO 35 J=1,I1T 

Zl eV (15d) 

CALDER PRA IAGIR, IT DEM C.K eV) 
V (NL, NL+IR) =V (NL, NL+IR) -MR 
V (NL+IR, NL) =V (NL+IR, NL) +MR 
V (NL+2*IR, NL) =V (NL+2*IR, NL) +MR*«DT 
V (NL+2*IR, NL+IR) =V (NL+2*IR, NL+IR) +0. 2*MR*xDT**2 
V (NL+2*IR, NL+2*IR) =V (NL+2*IR,NL+2* IR) +MR*DT**3/60. 


Determine vector of nodal acceleration at time t=0 
to enable starting procedure for reccurence scheme 


DO 50 I=1,IT 
Q(I)=0.0 
WVO (I) =0.0 
DO 60 I=1,IR 
WV1 (I) =0.0 
WV 1 (NL) =PYI 
IDGT=0 
IRHS=1 
CALL LEQT2F(M, IRHS,IR,IR,WV1,IDGT,WK, IER) 
DO 70 I=1,IR 
WVO (IT-IR+I1) =WV1 (I) 


Direct integration of a system equations of motion 
(using time-stepping scheme) and printing out the 
Fesuitsiy whichieerOmemtne sake- of Simplicity, /is 
limited only to displacement envelope at a midspan 
of a beam. Subroutine FORCE determines dynamic oil 
film force (vertical component, Fy) at each of IP 
time sub-intervals and subroutine SPLINT evaluates 
required integrals. 


NPT=0 
PREV=-0. 1 
LAST=0.0 
T=0.0D+00 
WRITE (6,3300)SMN,STP,UC,RSP,WT,PAR 
WRITE (6, 1109) 
DO 80 L=1,NDT 
IF(L .LE. NDT2) MAR=2 
IF(L .GT. NDT2) MAR=3 
CALL FORCE(T,DX,IP,TA,EPS,PSI,VEP,VPS,OMR,OMI,M 
WT,MC,ME,PAR,TX,PX,PY,OMT,OMS,OMA) 
CALL SPLINT GIP ,DXsTX, P¥sBB,CC,F1,F2,F3) 
Q(NL) =F 1 
Q(NL+IR) =F2 
Q(NL+2*IR) =F3 
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CALL IMNUEECVeWV0;, 1T; ITs 1, WV 1) 
BOWS 2=1 TT 
85 WV2 (1) =Q (1) -WV1 (I) 
CALL IMMULT (25 WV231LT TT #13 WVO) 
IF(LAST .GT. PREV .AND. LAST .GT. WVO(NL))GO TO 
+ 977 
GO TO 978 
977 TIME=T-DT 
NPT=NPT+ 1 
X=TIME/TO 
Y=LAST/YO 
WRITE (6,1110)NPT,TIME,LAST,X,Y 
978 PREV=LAST 
LAST=WVO (NL) 
80 CONTINUE 
NPT=999 
WRITE(6,1110)NPT 


1000 FORMAT (615) 

2000 FORMAT (1015) 

1001 FORMAT (1I6,5D20.8) 

1002 FORMAT (9D20.8) 

1100 FORMAT(1X,'Total number of finite elements used', 


+ Ss IS. 
te 1X, 'Number of finite elements in beam : 
i 8x, "=" pT3/ 
+ 1X, 'Number of finite elements in column ', 
+ 8X, '=',13/) 


1101 FORMAT (1X, 615) 
2101 FORMAT(1X,13,915) 


1102 FORMAT(1X,'Poisson ratio =',F19.5/ 

+ 1X,'Young modulus =',F19.5/ 

+ 1X, 'Density =',F19.5/ 

+ 1X, 'Shape factor ='F19.5// 

+ 1X, 'Total beam length ='F19.5/ 

+ 1X, 'Beam cross-section area =',F19.5/ 

+ 1X, 'Beam second moment of area =',F19.5/) 
1104 FORMAT(1X,'System 1st natural frequency=',F19.5/ 

+ 1X, 'System 1st natural period =',F19.5/ 

te 1X, 'System highest frequency =',F19.5/ 

+ 1X, 'System shortest period =',F19.5) 
1105 FORMAT (1X/16X,'STRUCTURE AND MATERIAL PROPERTIES ' / 

+ 16K, moore Up) 
1106 FORMAT (1X/1X, 'Array NG(NE,IE):'/ 

+ ee 578) 
1107 FORMAT (1X/16X,'FORCED VIBR. TRANSIENT SOLUTION'/ 

st. ot MTS Ae wh 
1108 FORMAT(1X,'Time step =',F19.5/ 

+ 1X,'No. of time steps may, 

+ 1X, 'Alpha =',F19.5/ 

+ 1X, 'Beta =',F19.5/ 

+ 1X, 'Gamma =',F19.5/ 
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+ 1kel Delta ='\F19.5/) 
1109 FORMAT (12X, 'TIME', 6X, 'DEFLECTION',8X,'TIME/TO', 

+ 9X, 'DRFactor'/ 

: (2k gh =o WPOX bee oe eh Stella s 

: 2E So 11) 
1110 FORMAT (1X,13,5X,F8.4,5X,E11.5,5X,F10.3,5X,F10.3) 
1111 FORMAT(1X,'Number of degree of freedom per', 

+ 1X, 'element mee y 

1X, -TOtaL*NO. Of DVOsFS~- of -constrained”; 

+ ixX-'system — =)e13/ 

a 1X, 'Order of finite time element', 

= 1X, 'matrices =v 13/ 

ae 1X, 'Nodal number corresponding to applied’, 

- ixXstload ='!,13/) 
1114 FORMAT(16X,'BEAM SIMPLY SUPPORTED AT ITS ENDS ' / 

fi Ley 52a ny) 
ads FORMAT (1X, 'COMPLEX TIMOSHENKO BEAM ELEMENT (TM544) 

ae se eA) 
1116 FORMAT(1X,'SIMPLE TIMOSHENKO BEAM ELEMENT' / 

ae be rh) 


1117 FORMAT (1X,16,3X,F15.8) 
1119 FORMAT (1X/) 


1122 FORMAT(1X,'Mass of the rotor =',F19.5/ 
oF 1X, 'Eccentricity of the mass =' F19.5/ 
+ 1X, 'Rotor/support mass ratio ='\F19.5/) 
1123 FORMAT(1X,'System ist natural frequency=',F19.5/ 
+ 1X,'System ist natural period =',F19.5/ 
a 1X, 'System highest frequency =' F19.5/ 
ae 1X, 'System shortest period =',F19.5/ 
+ 1X,'Time step/shortest period =',F19.5/ 
s 'O', 'Rotor operating speed =' F19.5/ 
+ 1X, 'Rotor acceleration time =',F19.5) 
1125 FORMAT(1X,'Rotor initial speed =' F19.5/) 
3300 FORMAT (5X, 'SMN=',F10.3,3X, 'STP=',F10.3,3X, 'e/C=',F 
+ 5X; 'NOy="Srevieok, 'W =',FS.1,5xX, 'PAR=",F10 
RETURN 
END 
ele 
ccc Subroutine FRAME (version for Model 2) 
ele 


SUBROUTINE FRAME(NBE,NCE,NE,IE,IR,IT,NL,ISM,MWK,M, 
- K,C,E,AA,BB,DD,WK,NG,Z,V,Q,WVO,WV1,WV2) 


REAL*8 M(IR,IR),K(IR,IR),C(IR,IR) ,BB(ISM) ,AA(ISM) 
REAL*8 DD(IR),E(IE,IE), Z(IT,IT),V(IT,IT) ,WK (MWK) 
REAL*8 WVO(CIT) ,WV1 (IT) ,WV2(IT) ,QCIT) 


REAL*8 ALP,BET,GAM,DEL,EPS,ZET,PI,GKA,BML,CNL,TO 
REAL¥S C1;C2,C3,C4,C5,CK.D1,D2,D3,D4,D5,F1,F2,F3 
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REAL*8 EL,EMO,SMO,ARE,SMA,RHO,PRO,ECC,MR,OMO,OMN 
REAL*8 OMR,OMT,OMS,OMA,DT,TA,T,TX,FX,PREV,LAST,X 
REAL*8 Y,X0,YO,PREVX,LASTX,PREVY,LASTY,HOM,TR,TN 
REAL*8 TG(4) ,HG(4) ,OMI,SHF 


INTEGER NG(NE, IE) 


Gaussian quadrature internal points and weighting 
factors 


TG (1) =-0.861136311594053D+00 
TG (2) =-0. 33998 1043584856D+00 
TG (3) =-TG (2) 

TG (4) =-TG (1) 

HG (1) =0.347854845137454D+00 
HG (2) =0.652145154862546D+00 
HG (3) =HG (2) 

HG (4) =HG (1) 

READ (5,1001)NDT,DT,ALP, BET, GAM, DEL 
IF(NDT .GT. 0) GO TO 449 
WRITE (6, 1112) 

GO TO 450 

WRITE (6, 1107) 

CONTINUE 

WRITE (6, 1113) 

IF(IE .EQ. 6) GO TO 453 
WRITE (6, 1115) 


GO TO 454 

WRITE (6, 1116) 

CONTINUE 

WRITE (6, 1100) NE, NBE, NCE 
WRITE (6,1111) IE,IR,IT,NL 


WRITE (6, 1106) 

TFALEWSEOUs 10), GOLTORSS A 

READ (5, 1000) ((NG(I,J) ,J=1, IE) , I=1,NE) 
WRITE (6,1101) ((NG(I,J) ,J=1, IE) , I=1,NE) 
GO TO 552 

READ (5,2000) ((NG(1,J) ,J=1, IE) , I=1,NE) 
WRITE (6,2101) ((NG(I,J) ,J=1, IE) ,1=1,NE) 
CONTINUE 

WRITE (6, 1105) 

READ (5, 1002) PRO, EMO,RHO,SHF,ECC,EPS 
PI=3.141592654D+00 

SMO=0.5*EMO/ (1.0+PRO) 

RHO=RHO/ 386.16 

READ (5, 1002) BML,ARE,SMA 

EL=BML/NBE 

READ (5, 1002)MR 

WRITE (6, 1102) PRO, EMO, RHO, SHF,BML,ARE, SMA 
GKA=SMO*SHF*ARE 
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Evaluate coefficients for static test 


C1=(BML**3) / (96.0*EMO*SMA) 
C2=(9.0*BML**«2) /32.0 
C3=BML/ (4.0*GKA) 
C4=0.5/EMO 
C5=(9.0*BML**3) / (64.0*EMO*ARE) 
CK=SMA/BML 
D1=1.0/(12.0*xEMO) 

D2=0.5 

D3=1.0/ (GKA*BML) 

D4=2.0/ (EMO*BML**2) 
D5=BML/ (4.0*EMO*ARE) 


Evaluate beam element mass and stiffness matrices 
and assemble them into system global matrices 


DO 10 I=1,IR 
DO 10 J=1,IR 

M(I,J)=0.0 

K(I,J)=0.0 
IF(IE .EQ. 6) GO TO 661 
CALL KCTMBM (EL, GKA,ARE,EMO,SMA,E) 
GO TO 662 
CALL KSTMBM (EL, GKA, ARE, EMO, SMA, E) 
CONTINUE 


DO 12 NK=1,NBE 
CALL ASSZW(E,NK,IE,NE,IR,K,NG) 
IF(IE .EQ. 6) GO TO 663 
CALL MCTMBM (EL, ARE, SMA, RHO,E) 
GO TO 664 
CALL MSTMBM (EL, GKA, ARE, EMO, SMA, RHO,E) 
CONTINUE 
DO 14 NK=1,NBE 
CALL ASSZW(E,NK,IE,NE,IR,M,NG) 


Read in column geometry parameters 

CNL - length of column 

ARE - cross-~secyional area 

SMA - second moment of cross-section area 


READ (5, 1002) CNL,ARE,SMA 
EL=CNL/NCE 
MR=MR+EPS*2.0*CNL*ARE*RHO 
GKA=SMO*SHF*ARE 


Continue evaluation of coefficients for analytical 
determination of maximum static displacements at 
midspan of a frame beam. 
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CK=CK*CNL/SMA 

C1=C1*(1.0+2.0%*CK) / (2.0+CK) 

C2=C2/ ( (CNL*GKA) * (2. O+CK) **2) 

C4=C4*CNL/ARE 

C5=C5/ (CNL* (2.0+CK) ) **2 
D1=D1*(2.0+3.0*CK) *CNL**3/ (SMA* (1.0+6.0%*CK) ) 
D2=D2*CNL/GKA 

D3=D3* (3.0*xCK*CNL/ (1.0+6.0%*CK) ) **2 

D4=D4* (CNL/ARE) * (3.O0*CK*CNL/ (1.0+6.0%*CK) ) **2 
P=1000.0 

YO=—-Px (C1+C2+C3+C4+C5) 

XO=P* (D1+D2+D3+D4+D5) 


Evaluate columns' elements matrices (stiffness and 
mass matrices and assemble them into system ones. 


IF(IE .EQ. 6) GO TO 666 

CALL KCTMCN (EL, GKA, ARE, EMO,SMA,E) 
GO TO 667 

CALL KSTMCN (EL,GKA,ARE,EMO,SMA,E) 
CONTINUE 


NS=NBE+ 1 

NF=NS+NCE- 1 

DO 16 NK=NS,NF 

CALL ASSZW(E,NK,IE,NE,IR,K,NG) 

IFC(IE .EQ. 6) GO TO 668 

DO 17 I=1,10 
BG?, D=-EG1, 1) 
E(6,1)=-E(6,1) 
E Ch, '-2 7, 1) 
E(1,6)=-E(1I,6) 
E(2,1)=-E(2,1) 
Ba, L)=-E (751) 
E(1,2)=-E(I,2) 
E(1,7)=-E(I,7) 

CONTINUE 

NS=NF+1 

NF=NS+NCE- 1 

DO 18 NK=NS,NF 

CALL ASSZW(E,NK,1E,NE,IR,K,NG) 

IF(IE .EQ. 6) GO TO 669 

CALL MCTMCN (EL, ARE, SMA, RHO,E) 

GO TO 670 

CALL MSTMCN (EL,GKA,ARE,EMO,SMA,RHO,E) 

CONTINUE 
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NS=NBE*+1 
NF=NS+NCE- 1 
DO 20 #£=°/\NK=NS,NF 
20 CALL ASSZW(E,NK,IE,NE,IR,M,NG) 
TRGiEe BO 26) “GO: "TO. 671 
DONZ4ra= 1, 10 
E(1,I)=-E(1,1) 
EG, ))=-8(6,1) 
BC 0) =-E (isi) 
E(1,6)=-E(1,6) 
E(2,1)=-E(2,1) 
BG, 0) == BiC7, 1) 
E(I,2)=-E(I, 2) 


ox E(1, 7) =-E(!, 7) 
671 CONTINUE 
NS=NF+1 
NF=NS+NCE- 1 


DO 22 NK=NS,NF 

22 CALL ASSZW(E,NK,IE,NE,IR,M,NG) 
LE (NDIG{CT<60) GO TO 772 
NX=NL-1 
M (NX, NX) =M (NX,NX)+MR 
WRITE(6, 1103) CNL,ARE,SMA 

771 CONTINUE 
M(NL,NL)=M(NL,NL)+MR 
EPS=MR/ (BML*ARE*RHO) 
WRITE (6, 1122)MR,ECC,EPS 


Determine system natural frequencies & eigenvector 


DO 37 I=1,IR 
JR=1 
DO 38 J=1,JR 
IS=I* (1-1) /2.0+0.5+d 
AA (IS) =K(1,J) 
BB(IS)=M(1,J) 
38 CONTINUE 
37 CONTINUE 
CALL EIGZS(AA,BB,IR,1,DD,C,IR,WK, IER) 
DO 39 I=1,IR 
39 DD (I) =DD(I) **0.5 
OMO=DD (1) 
HOM=DD (IR) 
TO=(PI*2.0D+00) /OMO 
TN=2.0*PI/HOM 
WRITE (6, 1104) OMO,TO,HOM, TN 
WRITE (6, 1120) 
DO 40 I=1,IR 
40 WRITE (6, 1117) 1,DD(I) 
DO 45 I=1,IR 
WRITE (6,1127)1I, (C(I,J) ,J=1,5) 
45 CONTINUE 
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FORMAT (1X,13, 10(2X,E10.4)) 
Static test 


IDGT=0 
DO 69 I=1,IR 
DO 69 J=1,IR 
CCLRJ)=K (1,3) 
DO 71 I=1,IR 
DD(I)=0.0 
DD (NL) =-P 
CALL LEQT2F(C,1,IR,IR,DD, IDGT,WK, IER) 
WRITE (6, 1119) 
DO 73 I=1,IR 
WRITE(6,1117)1,DD(1) 
ERR=DABS (100.0* (DD (NL) -YO) /YO) 
WRITE (6,1118)YO,ERR 
DO 74 I=1,IR 


DD(I)=0.0 
DD (NX) =1000.0 
IDGT=0 


DOLZSPTEILLIN 
DO 75 J=1,IR 
C(r;.J)=K (193) 
CALL LEQT2F(C,1,IR,IR,DD,IDGT,WK, IER) 
WRITE (6,1119) 
DO 78 I=1,IR 
WRITE(6,1117)1,DD(I) 
ERR=DABS (100.0* (DD (NX) -X0) /X0) 
WRITE(6,1121)X0,ERR 
CONTINUE 


READ (5, 1001) 1,OMN,OMO,HOM 

OMR=DEL*OMN 

YO= (MR*ECC*OMR**2) * (C1+C2+C3+C4+C5) 
TO=2*PI/OMO 

TN=2*PI/HOM 

TA=GAM*«TO 

ZET=DT/TN 

WRITE (6,1126) I 

WRITE (6, 1108) DT,NDT,ALP, BET, GAM, DEL 

WRITE (6,1124)OM0,TO,HOM,TN,ZET,OMN,OMR,TA 


Evaluate system global damping matrix 
DO 90 I=1,IR 
DO 90 J=1,IR 
C(I,J)=ALP*«M (I,J) +BET*K (I,J) 


Call subroutines to create submatrices of a finite 
time formulation and modify them. 
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CALESETZ12 (IR, ID. DT.MsC.K,Z) 
Z(NL,NL+IR)=Z(NL,NL+IR)+MR 
Z(NL+IR, NL) =Z(NL+IR, NL) -MR 

Z(NLt+IR,NL+IR) =Z(NL+IR,NL+IR) +MR*DT 
Z(NL+2*IR, NL) =Z(NL+2*IR, NL) -MR*DT 
Z(NL+2*IR,NL+IR) =Z (NL+2*IR,NL+IR) +O. 8*MR*DT*%*2 
Z(NL+2*IR,NL+2*IR) =Z (NL+2*IR, NL+2*IR) +MR*DT**3/60. 
Z (NX, NX+IR) =Z (NX, NX+IR)+MR 

Z (NX+IR,NX) =Z (NX+IR,NX)—-MR 

Z (NX+IR,NX+IR) =Z (NX+IR,NX+IR) +MR*DT 
Z(NX+2*IR, NX) =Z(NX+2*IR, NX) -MR*DT 

Z(NX+2*IR, NX+IR) =Z (NX+2*IR, NX+IR)+0.8*xMR*xDT**2 
Z(NX+2*IR, NX+2*IR) =Z (NX+2*IR, NX+2*1IR) +MR*DT**3/60. 
IDGT=0 


Invert modified submatrix (212) 


CALL LINV2F(Z,IT,IT,V,IDGT,WK, IER) 
DO 95 I=1,IT 

DO 95 J=1,I1T 

Z(1,J)=v(1,J) 

CALBSETZii CIR IT SDT. MC RV) 
V (NL, NL+IR)=V (NL,NL+IR) -MR 
V (NL+IR, NL) =V (NL+IR,NL)+MR 
V (NL+2*IR, NL) =V (NL+2*IR, NL) +MR*DT 
V (NL+2*IR, NL+IR) =V (NL+2*IR, NL+IR) +0. 2*MR*DT*x*2 
V (NL+2*IR, NL+2*IR) =V (NL+2*IR, NL+2*IR) +MR*DT**3/60. 
V (NX, NX+IR)=V (NX, NX+IR)-MR 
V (NX+IR, NX) =V (NX+IR,NX)+MR 
V (NX+2*IR, NX) =V (NX+2*IR, NX) +MR*DT 
V (NX+2*IR, NX+IR) =V (NX+2*IR, NX+IR) +O. 2*MR*DT**2 
V (NX+2*IR, NX+2*IR) =V (NX+2*IR, NX+2* IR) +MR*DT**3/60. 


Determine vector of nodal acceleration at time t=0 


DO 100 I=1,IT 
Q(1I)=0.0 
WVO (I) =0.0 
DO 110 I=1,IR 
WV1(I)=0.0 
WV1 (NL) =-MR*ECC*2.0*OMR/TA 
IDGT=0 
CALL LEQT2F(M,1,IR,IR,WV1,IDGT,WK, IER) 
POMIZ0 (IS1SIR 
WVO (IT-IR+I1) =WV1 (I) 


Direct integration of a system equations of motion 


NETA=O 
PREVA=—O7) 1 
LASTX=0.0 
NPeY=0 
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729 0259 PREVY=-0. 1 

730 0260 LASTY=0.0 

731 0261 T=0.0 

732 0262 WRITE (6, 1109) 

733 0263 DO 130 L=1,NDT 

734 0264 CALL LOADY(T,DT,TA,MR,OMR,ECC,TG,HG,F1,F2,F3) 
735 0265 Q (NL) =F 1 

736 0266 Q(NL+IR) =F2 

737 0267 Q(NL+2*IR) =F3 

738 0268 CALL LOADX(T,DT,TA,MR,OMR,ECC,TG,HG,F1,F2,F3) 
739 0269 Q (NX) =F 1 

740 0270 Q (NX+IR) =F2 

741 0271 Q(NX+2*IR) =F3 

742 O27Z CALL MMULT(V,WVO,IT,IT,1,WV1) 

743 0273 DO 135,1=791ft 

744 0274 135 WV2 (I) =Q (I) -WV1 (I) 

745 0275 CARL SMMULT (ZEWVZ SITS ITs 1 WVG) 

746 0276 IF(LASTY .GT. PREVY .AND. LASTY .GT. WVO(NL)) 
747 + GO TO 1977 
748 0277 GO TO 1966 

749 0278 1977 NPTY=NPTY+1 

750 0279 X=T/TO 

751 0280 Y=LASTY/YO 

752 0281 WRITE (6,1110)NPTY,T,LASTY,X,Y 

753 0282 1966 IF (LASTX .GT. PREVX .AND. LASTX .GT. WVO(NX)) 
754 “ GO TO 1988 
755 0283 GO TO 1955 

756 0284 1988 NPTX=NPTX+1 

757 0285 X=T/TO 

758 0286 Y=LASTX/YO 

759 0287 WRITE C7, 1110) NPTX, T;LASTX /X5Y 

760 0288 1955 T=T+DT 

761 0289 PREVX=LASTX 

762 0290 LASTX=WVO (NX) 

763 0291 PREVY=LASTY 

764 0292 LASTY=WVO (NL) 

765 0293 130 CONTINUE 

766 0294 NPTX=999 

767 0295 NPTY=999 

768 0296 WRITE (7, 1110) NPTX 

769 0297 WRITE (6, 1110) NPTY 

770 e 

Wht 0298 1000 FORMAT (615) 

WIZ 0299 2000 FORMAT (1015) 

sy 0300 1001 FORMAT (16,5D20.8) 

774 0301 1002 FORMAT (9D20.8) 

I5 0302 1100 FORMAT(1X,'Total number of finite elements used', 
776 + 8X, '=',13/ 

BI + 1X, 'Number of finite elements in beam a 
778 us ex, '=',13/ 

779 + 1X, 'Number of finite elements in column ', 
780 + ex, '=',13/) 
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1101 FORMAT (1X,615) 
2101 FORMAT(1X,13,915) 


1102 FORMAT(1X,'Poisson ratio =',F19.5/ 
3 1X,'Young modulus =',F19.5/ 
ts 1X, ‘Density ='F19.5/ 
te 1X, 'Shape factor =',F19.5// 
+ 1X, 'Total beam length =',F19.5/ 
+ 1X,'Beam cross-section area ='F19.5/ 
+ 1X,'Beam second moment of area =',F19.5/) 

1103 FORMAT (1X, 'Column height =',F19.5/ 
3 1X, 'Column cross-section area =',F19.5/ 
+ 1X,'Column second moment area =',F19.5/) 

1104 FORMAT(1X,'System ist natural frequency =',F19.5/ 
+ 1X, 'System 1st natural period =',F19.5/ 
ze 1X, 'System highest frequency Beate 
+ 1, ‘System shortest period ,F19.5) 

1105 FORMAT (1K/ 16K, "STRUCTURE AND MATERIAL PROPERTIES ' / 
- iC, Wigscac pease aseaucend co ae 18. 

1106 FORMAT (1X/ 1X, ‘Array NG(NE,IE):'/ 

+ eta '/) 

1107 FORMAT(1X/16X,'FORCED VIBR. TRANSIENT SOLUTION'/ 
+ See SX TE ye 

1108 FORMAT(1X,'Time step = ',F19.5/ 
a 1X,'No. of time steps eo 
+ 1X, 'Alpha = ' F19.5/ 
* 1X, 'Beta = ',F19.5/ 
+ 1X, 'Gamma = ',F19.5/ 
+ 1X, ‘Delta = ' | F19.5/) 

1109 FORMAT (12X, 'TIME',6X,'DEFLECTION' , 8x, 'TINE/TO", 

+ Gx: DRFactor' / 
+ 12 ee 6 Oke ee a 
2 OO Sais Gas a) 
1110 BORMAT (AX NIS¥5X SPSe4e5x Bi 45,5X,F10.3,5X,F 10. 3) 
1111 FORMAT(1X,'Number of degree of freedom per' 
+ 1X, ‘element pL ys 
a {X, “Total Noweof D.O.Fs of constrained’, 
+ ix)" system S='.13/ 
+ 1X, 'Order of finite time element', 
+ 1X, 'matrices =n / 
+ 1X, 'Nodal number corresponding to applied’, 
+ 1X2 LOaAdy= = 41S) 

1112 FORMAT (1X/16X, 'NATURAL FREQUENCIES & STATIC TEST'/ 
wu b+ bed SK! .) 

1113 FORMAT (16X,'RIGID PORTAL FRAME WITH CLAMPED LEGS'/ 
ae bt leap x! m 

1115 FORMAT(1X, ‘COMPLEX TIMOSHENKO BEAM ELEM. (1TM544) '/ 
ae eon wah) 

1116 FORMAT (1X, 'SIMPLE TIMOSHENKO BEAM ELEMENT’ / 

+ wt 7 fs) 


1117 FORMAT (1X,16,3X,F15.8) 
1118 FORMAT('O','static test (y-deflec. at midspan) '/ 
+ Ga uf 
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833 + 1X, 'exact solution =',F15.8/ 

834 + 1X, 'percentage error=',F15.8/) 

835 0321 1119 FORMAT (1X/) 

836 ©222 1120 FORMAT(1X/1X,'System natural frequencies:'/ 

837 a cae 47) 

838 0323 1121 FORMAT('0', 'static test (x-deflec. at midspan) '/ 
839 ue hwo! ih 

840 + 1X,'exact solution =',F15.8/ 

841 + 1X, 'percentage error=',F15.8/) 

842 0324 1122 FORMAT(1X,'Mass of the rotor = ' F19.5/ 
843 + 1X, ‘Eccentricity of the mass = ',F19.5/ 
844 + 1X, 'Rotor/support mass ratio = ' F19.5/) 
845 0325 1124 FORMAT(1X,'System 1st natural frequency= ',F19.5/ 
846 = 1X, 'System ist natural period = ',F19.5/ 
847 + 1X, 'System highest frequency = ' F19.5/ 
848 oF 1X, 'System shortest period =O! FA 25/ 
849 + 1X,'Time step/shortest period = ',F19.5/ 
850 e 'O', 'Highest passed frequency = ' F19.5/ 
851 e 1X, ‘Rotor operating speed = ',F19.5/ 
852 + 1X, 'Rotor acceleration time = ',F19.5/) 
853 0326 1125 FORMAT(1X,'Rotor initial speed = ',F19.5/) 
854 0327 1126 FORMAT(1X/1X,'ROTOR SPEED PASSING' ,I13,2X, 'NATUR', 
855 + 'AL FREQUENCY ') 

856 (e 

857 0328 RETURN 

858 0329 END 

859 é 

860 Ce 

861 ewe: Subroutine to multiply matrices 

862 Re 

863 C 

864 0001 SUBROUTINE MMULT(A,B,M,KK,N,C) 

865 0002 REAL*8 A(M,KK) ,B(KK,N) ,C(M,N) 

866 0003 DOs Lai iM 

867 0004 DO) li ena 

868 0005 C(1I,J)=0.0 

869 0006 DO 1 L=1,KK 

870 0007 1 Chg) =C Gd) +h GES) Bd) 

871 0008 RETURN 

SZ 0009 END 

Sts c 

874 ele: 

875 ASS Subroutine to create finite time element matrix 
876 Cee from system mass, damping and stiffness matrices 
877 ore (fone natrix. 212) 

878 G 

879 0001 SUBROUTINE: FRZt2 CIR ATT <DTeAyB, CZ) 

880 0002 REAL*8 A(IR,IR),BC(IR,IR),C(IR,IR),Z(IT,IT) ,DT 

881 0003 DO 1 I=1,IR 

882 0004 DO a et ky 

883 0005 Z(1,J)= Ei an). DTG, Dei 2 . 


884 0006 Z(I,J+IR)= BGO) oDTee 2x6 (1) 10. 
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Z(1,Jd+2*IR)= DUes3sxc (1 ,J)/ 120. 

Z(I+IR,J)= =A(T,J)+DTx(B(C1,J) /2.0tDT* 
Clr) the 14.) 

Z(I+IR,J+IR)= DE ACE sd) +DT* (2G) /10.— 


DECC, J) His. 7 oO.) 
WUsIR, J+2*IR)= DT**3*(-B(1,d)/120.+DT* 
CCE) / 2102) 
Z(I+2*IR,J)= DIX (GAC 53) +DT* (BG) * 
De) TeFDTXCCL J) 423. /84.)) 
Z(I+2*IR,JtIR)= (AC1,J)*4./5.+DTx(B(1,J) *13 
(105.-DTXC (1, 3) / 24.) ) *DTx*2 
Z(I+2*IR,J+2*IR) =(A(1I,J) /60.-DT* (B(I,J) /105. 
-DT*C (I,J) /336.) ) *DT**3 
CONTINUE 
RETURN 
END 


Subroutine to create finite time element matrix 
from system mass, damping and stiffness matrices 
(For, matrixes Zit) 


SUBROUTINE FTZ11(IR,IT,DT,A,B,C,V) 
REAE*GPAIK Re IR) sBCIR, PR) CCUIR, IR) VCIT, IT), DT 
DO 1 I=1,IR 

DOM tJ=1,IR 


V(I,J)= -B (UL U.EDIXC( Lt, 3) (2. 

V(1I,J+IR)= =K( ied) +DTxx2xC 1, J) / 10. 

V(1,J3+2*IR) = DTeA3*C CL) / 120. 

VCI+IR,J)= NO) tDT* Bi J) /2e+DT* 
C3) /7.) 

V (I+IR,J+IR)= DEA*2* (—B (10) /10.4+DT* 


CCie) 46/005.) 
VCItIR,J+2*IR)= DT**3*(-B(I,J)/120.+DTx 
Cie) 2605) 
V(I+2*IR,J)= DT#({A (1,3) +DT* (BCT, J) 
2./Te+DT*C (1, 0)%5./84.)) 
V(IF2xIR.0+IR)= SAGs) /52-DT* (BIg a) . / 
$05. -DT*C.(T d) /56.) )*DPx*2 
V (1+2*IR, J+2*IR) =(A(1,J) /60.-DT* (B(I,J) /140. 
-pT*C (I,J) /560.)) *DT**3 
CONTINUE 
RETURN 
END 


Subroutine to create element stiffness matrix for 
the simple Timoshenko beam element 


SUBROUTINE KSTMBM (EL, GKA, ARE, EMO, SMA, EK) 
REAL*8 EK(6,6) 
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938 
339 
940 
941 
942 
943 
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951 
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955 
954 
955 
956 
Soy / 
958 
goo 
960 
961 
962 
963 
964 
965 
966 
367 
968 
969 
S10 
OG 
OZ 
ATS 
974 
O75 
976 
QT 
978 
Si 
980 
98 1 
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987 
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REAL*8 EL,GKA,ARE,SMA,EMO,FI,C 
FI=12.*EMO*SMA/ (GKA*EL**2) 
C= (EMO*SMA/EL**3) /(1.0+FI) 
EK (1, 1) =ARE*xEMO/EL 
EK(2, 1) Oro 
EK (2,2) =12.*C 
EK (3,1)=0.0 
EK (3,2) =6.0*xEL*C 
EK (3,3) =C* (4. + FI) *EL*«*2 
EK (4,1) =-EK (1,1) 
EK (4,2)=0.0 
EK (4,3) =0.0 
EK (4, 4) =EK (1, 1) 
EK (5,1)=0.0 
EK (5, 2) =-EK (2, 2) 
EK (5,3) =-EK (3,2) 
EK (5,4) =0.0 
EK (5,5) =EK (2, 2) 
EK (6, 1)=0.0 
EK (6, 2) =EK (3, 2) 
EK (6, 3) =C* (2.0-FI) *EL**2 
EK (6,4) =0.0 
EK (6, 5) =EK (5, 3) 
EK (6, 6) =EK (3, 3) 
DO 10 I=1,5 

K=I+1 

DO 10 J=K,6 

EK (I,J) =EK (J, 1) 

RETURN 
END 


Subroutine to create element mass matrix for the 
Simple Timoshenko beam element 


SUBROUTINE MSTMBM (EL, GKA, ARE, EMO, SMA, RHO, EM) 

REAL*8 EM(6,6) ,EL,GKA,ARE,EMO,SMA,RHO,FI,C,D 

FI=12.*EMO*SMA/ (GKA*EL**2) 

C= (RHO*ARE*EL) / (1. O+FI) **2 

D= (SMA/ARE) /EL**2 

EM (1, 1) =RHO*ARE*EL/3.0 

EM(2,1)=0.0 

EM (2,2) =C*(13./35.+0. 7*FI+FI**2/3.+1.2*D) 

EM(3,1)=0.0 

EM (3, 2) =C*xBL* (11.5/210.5411.*F1/120.+F1**2/24.+D* 
(0.1-FI/2.)) 

EM (3,3) =C*xEL**2*(1./105.+FI/60.+FI**2/120.+D* 
(20/15.0F1/6.781e42/3.0)) 

EM(4,1)=EM(1,1)/2. 

EM (4,2)=0.0 

EM (4,3) =0.0 
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1000 
1001 
1002 
1003 
1004 
1005 
1006 
1007 
1008 
1009 
1010 
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1016 
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1020 
1021 
1022 
1023 
1024 
1025 
1026 
1027 
1028 
1029 
1030 
1031 
1032 
1033 
1034 
1035 
1036 
1037 
1038 
1039 
1040 


0015 
0016 
0017 
0018 


0019 
0020 
0021 
0022 
0023 


0024 
0025 
0026 


0027 
0028 
0029 
0030 
0031 
0032 


0001 
0002 
0003 
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0006 
0007 
0008 
0009 
0010 
0011 
0012 
0013 
0014 
0015 
0016 
0017 
0018 
0019 
0020 
0021 
0022 
0023 
0024 
0025 


ccc 


+ 


EM (4,4) =EM (1, 1) 
EM(5,1)=0.0 


EM (5,2) =C*«(9./70.+0. 3*FI+F1I**2/6.-1.2*D) 
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EM (5,3) =C*EL* (13. /420.+3.*FI/40.+FI**2/24.-D* (0. 1 


-FI/2.)) 
EM (5,4) =0.0 
EM (5,5) =EM (2,2) 
EM (6,1) =0.0 


EM (6,2) =-EM (5, 3) 


EM (6, 3) =C*xEL**2% (Dx (FI*¥*2/6.-F1/6.-1./30.)-1./140. 


-FI/60.-F1I**2/120.) 


EM (6,4) =0.0 
EM (6,5) =-EM (3,2) 
EM (6,6) =EM (3, 3) 


DO 10 1=1,5 
K=I+1 
DO 10 J=K,6 
EM (1,J)=EM(J,1) 
RETURN 
END 


Subroutine to create element stiffness 
the simple Timoshenko beam element (for column) 


SUBROUTINE KSTMCN (EL, GKA, ARE, EMO, SMA, EK) 


REAL*8 EK(6,6) 


REAL*8 EL,GKA,ARE,SMA,EMO,FI,C 


FI=12.*EMO*xSMA/ (GKA*EL**2) 
C= (EMO*SMA/EL**3) /(1.0+FI) 
EK(1,1)=12.%*C 

EK (2,1) =0.0 

EK (2, 2) =ARE*EMO/EL 

EK (3, 1) =-6.0*EL*C 

EK (3,2) =0.0 

EK (3, 3) =C* (4. +FI) *EL**2 
EK (4, 1) =-EK(1, 1) 

EK (4,2) =0.0 

EK (4, 3) =-EK (3, 1) 

EK (4,4) =EK (1,1) 
EK(5,1)=0.0 

EK (5,2) =-EK (2,2) 
EK(5,3)=0.0 

EK (5,4) =0.0 

EK (5,5) =EK (2, 2) 

EK (6, 1) =EK (3, 1) 

EK (6,2) =0.0 

EK (6, 3) =C* (2.0-FI) *EL**2 
EK (6, 4) =-EK (3,1) 

EK (6,5) =0.0 
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1063 
1064 
1065 
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EK (6,6 

DO 10 
K=I 
DO 


RETURN 
END 
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simple 


SUBROU 
REAL*8 
FI=12. 
C= (RHO 
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) =EK (3, 3) 
dee a ae 

a | 

Ord =k, 6 


Bild =EK (J, 1) 


tine to create element mass matrix for the 
Timoshenko beam element (for column) 


TINE MSTMCN (EL, GKA, ARE, EMO, SMA, RHO, EM) 
EM (6,6) ,EL,GKA,ARE,EMO,SMA,RHO,FI,C,D 
*EMO*SMA/ (GKA*EL**2) 
*ARE*EL) / (1. O+FI) **2 


D= (SMA/ARE) /EL**2 


EM (1,1 
EM (2, 1 
EM (2,2 
EM (3, 1 


EM (3,2 
EM (3, 3 


EM (4, 1 
EM (4,2 
EM (4, 3 


EM (4,4 
EM (5, 1 
EM (5,2 
EM (5,3 
EM (5,4 
EM (5,5 
EM (6, 1 
EM (6,2 
EM (6,3 


EM (6,4 
EM (6,5 
EM (6,6 
DO 10 
K=I 
DO 


RETURN 
END 


)=0*(13./35.+0. 7XFI+FI**2/3.4+1.2*D) 

)=0.0 

) =RHO*xARE*EL/3.0 

)=-C¥EL* (11./210.+11.*FI/120.+F1I**2/24.+D* 
COMI-FI/ 2.) 

)=0.0 

J=CeEL**2* (1./105.+F1I/60.+F1I**2/120.+D* 
(22715 FP 1/6. + 
FI**2/3.0)) 

\=cx (9./70.+0.3*FI+FI**2/6.-1.2*D) 

)=0.0 

)=-C*¥EL* (13. /420.4+3.*F1I/40.+FI**2/24.—-Dx 
COc= F172) 

)=EM (1, 1) 

)=0.0 

)=EM (2,2) /2.0 

)=0.0 

)=0.0 

) =EM (2, 2) 

) =-EM (4, 3) 

)=0.0 

)=CHEL*4A2* (D* (FI*¥*2/6.-F1/6.-1./30.)-1./140. 
-F1/60.—-FI**2/120.) 

) =-EM (3, 1) 

)=0.0 

) =EM (3, 3) 

I=1,5 

+1 

10 J=K,6 

EM (I,J) =EM (J,I) 
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GEG Subroutine to calculate components (in X-direc.) 
CCC of load vector. Subroutine uses 4-point Gaussian 
CEE quadrature to evaluate integrals. 
ce 
Cc 
SUBROUTINE LOADX(T,DT,T1,MR,OMR,ECC,TG,HG, 
+ PISR2 253) 
REAL*8 T,DT,T1,MR,OMR,ECC,F1,F2,F3,TG(4) ,HG(4) , TX 
REAL*8 OMT,OMS,OMA, FX 
F1=0.0 
F2=0.0 
F3=0.0 
DO 10 I=1,4 
TX=DT* (TG (I) +1.0)/2.0 
PCr eer. T1) GOeTO 4 
OMT=OMR*T 1 * ((T+TX) /T1) **2* (3.0- (T+TX) /T1) /3.0 
OMS=OMR* (T+TX) * (2.0- (T+TX) /T1) /T1 
OMA=OMR*2.0* (1.0-(T+TX) /T1) /T1 
GO TO 2 
1 OMT=OMR* ( (T+TX) -T1/3.0) 
OMS=OMR 
OMA=0.0 
Z CONTINUE 
FX=MR*ECC* (OMS**2*DCOS (OMT) +OMA*DSIN (OMT) ) 
F1=F1+HG (1) *FX 
F2=F2+HG (1) *TX*FX 
10 F3=F3+HG (1) *TX**2*FX 
F1=F1*DT/2.0 
F2=F2*DT/2.0 
F3=F3*DT/2.0 
RETURN 
END 


Cee Subroutine to calculate components (in Y-direc.) 
CCE of load vector. Subroutine uses 4-point Gaussian 
Cee quadrature to evaluate integrals 


ce 
Cc 
SUBROUTINE LOADY(T,DT,T1,MR,OMR,ECC,TG,HG,F1,F2 
ae ,ES) 
fe: 
REAL*8 T,DT,T1,MR,OMR,ECC,F1,F2,F3,TG(4) ,HG(4) , TX 
REAL*8 OMT,OMS,OMA,FX 
C 
F1=0.0 
F2=0.0 
F3=0.0 
DO 10 I=1,4 


b) 
TX=DT* (TG(I)+1.0)/2.0 
TEOT+~GTs . TiRGOsTOM 
OMT=OMR*T 1% ((T+TX) /T1) **2%* (3.0- (T+TX) /T1) /3.0 


ap eV Ca tans) -0 ion 
hae Mee 
tits 


i 


Ro 
Wes 


Cea ae  eewmgod-tene) 2 


1145 
1146 
1147 
1148 
1149 
1150 
tid 
trS2 
TESS 
1154 
WSS 
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tIG3 
tT62 
1163 
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1165 
1166 
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hieg 
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1174 
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1178 
TTS 
1180 
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1183 
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tHES 
1186 
1187 
1188 
1189 
T1SG 
1191 
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bags 
1194 
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OMS=OMR* (T+TX) * (2.0- (T+TX) /T1) /T1 
OMA=OMR*2.0*(1.0-(T+TX) /T1) /T1 
GOFFORZ 
OMT=OMR« ( (T+TX) -T1/3.0) 
OMS=OMR 
OMA=0.0 
CONTINUE 
FX=MR*ECC* (OMS**2*DSIN (OMT) -OMA*DCOS (OMT) ) 
F1=F1+HG (1) «FX 
F2=F2+HG (1) *TX*FX 
F3=F3+HG (1) «TX**2*FX 
F1=F1*DT/2.0 
F2=F2*DT/2.0 
F3=F3*DT/2.0 
RETURN 
END 


Subroutine to create Timoshenko beam (beam-truss) 
element (TM544) stiffness matrix (for beam) 


SUBROUTINE KCTMBM (EL, GKA,ARE,EMO,SMA,K) 
REAL*8 K(10,10),EL,GKA, ARE, SMA, EMO 
K(1, 1) =ARE*EMO/EL 

K(1,2)=0.0 

K(1,3)=0.0 

K(1,4)=0.0 

K(1,5)=0.0 

K(1,6)=-K(1, 1) 

K(1,7)=0.0 

K(1,8)=0.0 

K(1,9)=0.0 

K(1,10)=0.0 

K (2,2) =6. *GKA/ (5.0*EL) 

K (23) =K (2, 2) *EL/2.0 

K (2,4) =-K(2,3) / (6. O*ARE) 
K(2,5)=-K (2,4) *ARE*EL/SMA 

K(2,6)=0.0 

K(257) =k (292) 

K (2,8) =K (2,3) 

K (2,9) =K(2,4) 

K(2, 10) =-K (2,5) 

K (3,3) =K (2,3) *EL/2.0+6.0*EMO*xSMA/ (5.0*EL) 
K (3,4) =K (2,4) *EL/2.0 

K(3,5)=K(2,5) *EL/2.0+EMO/10.0 
K(3,6)=0.0 

RG, 7) ==K(2,3) 

K (3,8) =K (2,3) *EL/2.0-6.0*EMO*SMA/ (5.0*EL) 
K (3,9) =K (3, 4) 
K(3,10)=-K (2,5) *EL/2.0+EMO/10.0 

K (4,4) =-4.0*K (3,4) *2.0/ (3.0*ARE) 
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1198 
+299 
1200 
1201 
1202 
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1204 
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1207 
1208 
1209 
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t2x3 
1214 
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1216 
1207 
1218 
t219 
1220 
f221 
t222 
t223 
1224 
1225 
1226 
1227 
1228 
229 
1230 
1231 
1232 
1233 
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K (4,5) =K (3,4) xEL/ (6.0*SMA) 
K(4,6)=0.0 
K (4,7) =-K (2,4) 
K (4,8) =K (3,4) 
K(4,9)=-K (4,4) /4.0 
K (4, 10) =-K (4,5) 
K(5,5)=K(2,5) *EL**2/(12.0*SMA)+2.*EMO*EL/ (15.*SMA) 
K(5,6)=0.0 
K (5,7) =-K(2,5) 
K (5,8) =-K(3, 10) 
K(5,9)=K(4,5) 
K(5,10)=-K (2,5) *EL**2/ (12.0*SMA) -EMO*EL/ (30.0*SMA) 
K(6,6)=K(1, 1) 
K(6,7)=0.0 
K (6,8) =0.0 
K(6,9)=0.0 
K(6,10)=0.0 
ECh, 7) =KC2;2) 
K(7,8) =-K(2, 3) 
K(7,9)=K (4,7) 
K(7,10)=K(2,5) 
K(8,8) =K (3, 3) 
K (8,9) =K (3, 4) 
K(8,10)=-K(3,5) 
K(9,9)=K (4,4) 
K (9, 10) =K (4, 10) 
K(10, 10) =K (5,5) 
DOS1O1T=2410 

N=I-1 

DO 10 J=1,N 

KGI $d) =K (WI, 1) 

RETURN 
END 


Subroutine to create Timoshenko beam (beam-truss) 
element (TM544) stiffness matrix (for column) 


SUBROUTINE KCTMCN (EL, GKA, ARE, EMO, SMA,K) 
REAL*8 K(10,10) ,EL,GKA,ARE,SMA, EMO 
K(1,1)=1.2*GKA/EL 

K(1,2)=0.0 

K(1,3)=-0.6*GKA 

K(1,4)=0. 1*GKA/ARE 

K(1,5)=-0. 1*GKA*EL/SMA 
K(1,6)=-K(1, 1) 

Kit, 7) =0.20 

K(1,8)=K(1,3) 

K(1,9)=K(1,4) 

K(1,10)=-K(1,5) 

K (2,2) =ARE*EMO/EL 
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1265 
1266 
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tZaa 
1272 
i273 
1274 
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1276 
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WASh| 
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1284 
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1287 
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1289 
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1294 
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0014 
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cc 


K (2,3) 
K(2,4) 
K(2,5) 
K (2,6) 
K (2,7) 
K(2,8) 
K(2,9) 
KC2,00 
K(3,3) 
K(3,4) 
K(3,5) 
K (3,6) 
K(3,'7) 
K (3,8) 
K (3,9) 
K(3,, 10 
K (4,4) 
K (4,5) 
K (4,6) 
K (4,7) 
K (4,8) 
K (4,9) 
K(4,10 


=0.0 

=0.0 

=0.0 

=0.0 

=-K (2,2) 

=0.0 

=0.0 

)=0.0 

=-K(1,3) *EL/2.0+6.0*EMO*xSMA/ (5.0*EL) 
=-K(1,4) *EL/2.0 

=-K(1,5) *EL/2.0+EMO/10.0 
=-K(1,3) 

=0.0 

=K (3,6) *EL/2.0-6.0*EMO*SMA/ (5.0*EL) 
=K (3,4) 

)=K (1,5) *EL/2.0+EMO/10.0 
=-4,0*K (3,4) *2.0/ (3.0*ARE) 
=K (3,4) *EL/ (6.0*SMA) 
=-K(1,4) 

=0.0 

=K (3,4) 

=-K (4,4) /4.0 

)=-K(4,5) 


153 


K(5,5)=-K (1,5) *EL**2/(12.0*SMA)+.4*EMO*xEL/ (3.*SMA) 


K(5,6) 
K(5,7) 
K(5,8) 
K(5,9) 


=-K(1,5) 
=0.0 

=-K (3,10) 
=K (4,5) 


K(5,10)=K(1,5) *EL**2/(12.0*SMA) -EMO*EL/ (30.0*SMA) 


K(6,6) 
K(6,9) 
K(6,8) 
K (6,9) 
K (6,10 
Kg 7) 
K(7,8) 
K(7,9) 
KG. 0 
K(8,8) 
K(8,9) 
K(8,10 
K(9,9) 
K(9,10 
KC108 4 
DO 10 
N=I 
DO 


RETURN 
END 


=K(1,1) 
=0.0 
=-K(1,3) 
=-K(1,4) 
)=K(1,5) 
=K (2,2) 
=0.0 
=0.0 
)=0.0 

=K (3,3) 
=K (3,4) 
)=-K (3,5) 
=K (4,4) 
)=K (4, 10) 
0) =K (5,5) 
I=2,10 

-1 

10 J=1,N 
K(1,J)=K(J,1I) 
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1332 
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ccc 


Subroutine to create Timoshenko beam element 
(TM544) consistent mass matrix (for beam) 


SUBROUTINE MCTMBM (EL, ARE,SMA,RHO,M) 
REAL*8 M(10,10) ,EL,ARE,SMA, RHO 
M(1,1)=RHO*ARE*EL/3.0 

M(1,2)=0.0 

M(1,3)=0.0 

M(1,4)=0.0 

M(1,5)=0.0 

M(1,6)=M(1,1)/2.0 

M(1,7)=0.0 

M(1,8)=0.0 

M(1,9)=0.0 

M(1,10)=0.0 

M (2,2) =RHO*ARE*EL* 13/35.0 

M (2,3) =RHO*xARE*EL**2* 17/280.0 
M(2,4)=-M (2,3) *44/ (51*ARE) 
M(2,5)=M (2,3) *EL*7/ (102*SMA) 
M(2,6)=0.0 

M(2,7)=M(2,2) *9/26.0 
M(2,8)=-M (2,3) *11/17.0 

M(2,9)=-M(2,4) *13/22.0 

M(2,10)=M(2,5) 

M(3,3)=M (2,3) *EL/4.5+13*RHO*SMA*EL/35. 
M(3,4)=M(2,4) *19*EL/88.0 

M(3,5)=M (2,5) *EL*11/42.+11*RHO*EL**2/210. 
M(3,6)=0.0 

M(3,7)=-M (2,8) 

M(3,8)=-M (2,3) *EL*28/153.+9*RHO*SMA*EL/70. 
M(3,9)=-M (3,4) *45/57.0 

M(3, 10) =M (2,5) *EL*11/42.-13*RHO*EL**2/420. 
M(4,4)=-M (3,4) *16/ (19*ARE) 
M(4,5)=-M(2,5) *9*EL/ (42*ARE) 
M(4,6)=0.0 

M(4,7)=M (2,4) *13/22.0 

M(4,8)=M (3,9) 

M(4,9)=-M (4,4) *3/4.0 

M(4,10)=M (4,5) 

M(5,5)=-M (4,5) * (EL*ARE/ (9*SMA) +32/ (3*EL) ) 
M(5,6)=0.0 

M(5,7)=M(2,5) 

M(5,8)=-M(3, 10) 

M(5,9)=-M (4, 10) 
M(5,10)=M(5,9) * (EL*ARE/ (9*SMA) -8/EL) 
M(6,6)=M(1,1) 

M(6,7)=0.0 

M(6,8)=0.0 

M(6,9)=0.0 

M(6,10)=0.0 

M(7,7)=M (2,2) 
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M(7,8)=-M(2,3) 
M(7,9)=-M (2,4) 
M(7,10)=M (2,5) 
M(8,8)=M (3,3) 
M(8,9) =M (3,4) 
M(8,10)=-M (3,5) 
M(9,9)=M(4,4) 
M(9,10)=M(5,9) 
M(10,10)=M(5,5) 
BO. S10 P=2,010 

N=I-1 

DO 10 J=1,N 

M(I,J)=M(J,1) 

RETURN 
END 


Subroutine to create Timoshenko beam element 
(T™544) consistent mass matrix (for column) 


SUBROUTINE MCTMCN (EL, ARE, SMA,RHO,M) 
REAL*8 M(10,10),EL,ARE,SMA,RHO 

M(1, 1) =RHO*ARE*EL*13/35.0 

M(1,2)=0.0 

M(1,3) =-RHO*ARE*EL**2*17/280.0 
M(1,4)=-M (1,3) *44/ (51*ARE) 
M(1,5)=M(1,3) *EL*7/(102*SMA) 
M(1,6)=M(1,1)*9/26.0 

M(1,7)=0.0 

M(1,8)=-M (1,3) *11/17.0 

M(1,9)=-M(1,4) *13/22.0 

M(1,10)=M(1,5) 

M(2,2) =RHO*xARE*EL/3.0 

M(2,3)=0.0 

M(2,4)=0.0 

M(2,5)=0.0 

M(2,6)=0.0 

M(2,7)=M(2,2)/2.0 

M(2,8)=0.0 

M(2,9)=0.0 

M(2,10)=0.0 

M(3,3)=-M(1,3) *EL/4.5+13*RHO*xSMA*EL/35. 
M(3,4)=-M (1,4) *19*EL/88.0 

M(3,5)=-M (1,5) *EL*11/42.+11*RHO*EL**2/210. 
M(3,6)=-M(1,8) 

M(3,7)=0.0 

M(3,8)=M (1,3) *EL*28/153.+9*RHO*SMA*EL/70. 
M(3,9)=-M (3,4) *45/57.0 

M(3,10)=-M (1,5) *EL*11/42.-13*RHO*EL**2/420. 
M(4,4)=-M (3,4) *16/ (19*ARE) 
M(4,5)=M (1,5) *9*EL/ (42*ARE) 
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20 
15 
10 


M (4,6) 
M (4,7) 
M(4,8) 
M(4,9) 
M(4,10 
M(5,5) 
M(5,6) 
M(5,7) 
M(5,8) 
M(5,9) 
M(5,10 
M (6,6) 
M (6,7) 
M (6,8) 
M(6,9) 
M(6,10 
WCU?) 
M(7,8) 
M(7,9) 
MCT ato 
M(8,8) 
M(8,9) 
M(8,10 
M(9,9) 
M(9,10 
MC1O 8 
DO 10 
N=I 
DO 


RETURN 
END 


Subrou 
global 


SUBROU 
REAL*8 
INTEGE 
DO 10 

DO 


RETURN 
END 


=M(1,4) *13/22.0 
=0.0 

=M (3,9) 

=-M (4,4) *3/4.0 
)=M (4,5) 

=-M (4,5) « (EL*ARE/ (9*SMA) +32/ (3*EL,) ) 
=M(1,5) 

=0.0 

=-M (3, 10) 
=-M (4, 10) 

)=M (5,9) * (EL*ARE/ (9*SMA) -8/EL) 
=M(1,1) 

=0.0 

=-M(1,3) 

=-M (1,4) 
)=N-Gr, 5) 

=M (2,2) 

=0.0 

=0.0 

)=0.0 

=M (3,3) 

=M (3,4) 

)=-M (3,5) 

=M (4,4) 

)=M (5,9) 

0) =M (5,5) 
I=2,10 

-1 

10 J=1,N 
M(I,J)=M(J,1) 


tine to assemble element matrices into 
ones 


TINE ASSZW(EMAT,NK,IE,NE,IR,GMAT,NSORT) 
EMAT (IE, IE) ,GMAT (IR, IR) 

R NSORT (NE, IE) 

be 4e 

15 J=1,IE 

N1=NSORT (NK, I) 

N2=NSORT (NK, J) 

EE CONG /SGT ER) SOR. (N2 .GT. IR))).GO TO-20 
GMAT (N1,N2) =GMAT(N1,N2)+EMAT (I,J) 

CONTINUE 

CONTINUE 

CONTINUE 


igh eee ri ae 
AB ee ONE 2 


1457 
1458 
1459 
1460 
1461 
1462 
1463 
1464 
1465 
1466 
1467 
1468 
1469 
1470 
1471 
1472 
1473 
1474 
1475 
1476 
1477 
1478 
1479 
1480 
148 1 
1482 
1483 
1484 
1485 
1486 
1487 
1488 
1489 
1490 
1491 
1492 
1493 
1494 
1495 
1496 
1497 
1498 
1499 
1500 
1501 
1502 
1503 
1504 
tSOS5 
1506 
1507 
1508 


0001 


0002 
0003 
0004 


0005 
0006 
0007 
0008 
0009 
0010 
0011 
0012 
0013 
0014 
0015 
0016 


0017 
0018 
0019 
0020 
0021 
0022 
0023 
0024 
0025 
0026 
0027 
0028 


0029 
0030 
0031 
0032 
0033 
0034 
0035 
0036 
0037 
0038 


155 


Subroutine FORCE to calculate components (in X and 
andeye directions)» of ‘dynamic oil-film forces in 
in journal bearing. S/r uses 4th order Runge-Kutta 
method to solve equations of journal centre motion 


SUBROUTINE FORCE(T,DX,IP,TA,EPS,PSI,VEP,VPS,OMR, 
OMI,MAR,WT,MC,ME,PAR,TX,PX,PY,OMT,OMS, OMA) 


REAL*8 T,DX,TA,EPS,PSI,VEP,VPS,OMR,OMI,OMT,OMS,OMA 
REAL*8 TX(IP) ,PX(IP) ,PY(IP) ,WIT,MC,ME,P1,P2 
REAL*8 AK,AL,AM,AN,GK,GL,GM,GN,S1,S2,S3,S4 


CALL ROTTSA(MAR,T,OMI,OMR,TA,OMT,OMS,OMA) 
IPM1=IP-1 
DO 10 J=1,IPM1 
CALL BF(PAR,EPS,VEP,VPS,OMS,P1,P2) 
TX (J) =T 
PX (J) =- (P1*DSIN (PSI) +P2*DCOS (PSI) ) 
PY (J) =P1*DCOS (PSI) -P2*DSIN (PSI) +WT 
S1=EPS 
S2=PSI 
S3=VEP 
S4=VPS 
CALL RKA(DX,P1,P2,WT,ME,MC,EPS,PSI,VEP,VPS,OMT, 
OMS ,OMA,AK,AL,AM,AN) 
GK=AK 
GL=AL 
GM=AM 
GN=AN 
EPS=S1+AM/2.0D+00 
PSI=S2+AN/2.0D+00 
VEP=S3+AK/2.0D+00 
VPS=S4+AL/2.0D+00 
T=T+DX/2.0D+00 
CALL ROTTSA (MAR,T,OMI,OMR,TA,OMT,OMS,OMA) 
CALL BF(PAR,EPS,VEP,VPS,OMS,P1,P2) 
CALL RKA(DX,P1,P2,WT,ME,MC,EPS,PSI,VEP,VPS,OMT, 
OMS ,OMA,AK,AL,AM, AN) 
GK=GK+AK*2.OD+00 
GL=GL+AL*2.0D+00 
GM=GM+AM*2.0D+00 
GN=GN+AN*2.0D+00 
EPS=S1+AM/2.0D+00 
PSI=S2+AN/2.0D+00 
VEP=S3+AK/2.0D+00 
VPS=S4+AL/2.0D+00 
CALL BF(PAR,EPS,VEP,VPS,OMS,P1,P2) 
CALL RKA(DX,P1,P2,WT,ME,MC,EPS,PSI,VEP,VPS,OMT, 
OMS ,OMA,AK,AL,AM,AN) 
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GK=GK+AK*2.0D+00 
GL=GL+AL*2.0D+00 
GM=GM+AM*2 .O0D+00 
GN=GN+AN*2.0D+00 
EPS=S1+AM 
PSI=S2+AN 
VEP=S3+AK 
VPS=S4+AL 
T=T+DX/2.0D+00 
CALL ROTTSA(MAR,T,OMI,OMR,TA,OMT,OMS,OMA) 
CALL BF(PAR,EPS,VEP,VPS,OMS,P1,P2) 
CALL RKA(DX,P1,P2,WT,ME,MC,EPS,PSI,VEP,VPS,OMT, 
OMS ,OMA,AK,AL,AM, AN) 
EPS=S 1+ (GM+AM) /6.0D+00 
PSI=S2+ (GN+AN) /6.O0D+00 
VEP=S3+ (GK+AK) /6.0D+00 
VPS=S4+ (GL+AL) /6.0D+00 
CONTINUE 
CALL BF(PAR,EPS,VEP,VPS,OMS,P1,P2) 
J=IP 
TX (J) =T 
PX (J) =- (P1*DSIN (PSI) +P2*DCOS (PSI) ) 
PY (J) =P1*DCOS (PSI) -P2*DSIN (PSI) +WT 
RETURN 
END 


Subroutine RKA (Runge-Kutta 4th order) to calcula- 
te increments in displacements and velocities at 
each time step 


SUBROUTINE RKA(DX,P1,P2,WT,ME,MC,EPS,PSI,VEP,VPS, 
OMT,OMS,OMA,AK,AL,AM, AN) 


REAL*8 DX,P1,P2,WT,ME,MC,EPS,PSI,VEP,VPS,OMT,OMS 
REAL*8 OMA,A,SP,CP,SA,CA,TW,AK,AL,AM,AN 


A=OMT-PSI 

SP=DSIN (PS1) 

CP=DCOS (PSI) 

SA=DSIN (A) 

CA=DCOS (A) 

TW=2.0D+00 

AK=DX* ((P1+WT*CP+ME*xOMS *«*2*CA+ME*OMA*SA) /MC+tEPS* 
VPS**2) 

AL=DX* ( (P2-WT*SP+ME*xOMS**2*SA-ME*OMA*CA) /MC-VEP* 
VPS*«TW) /EPS 

AM=DX*VEP 

AN=DX*VPS 

RETURN 

END 
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Subroutine BF to calculate bearing dynamic forces 
(in polar coordinates) at each time subinterval 


SUBROUTINE BF(PAR,EPS,VEP,VPS,OMS,P1,P2) 
REAL*8 PAR,EPS,VEP,VPS,OMS,P1,P2 
REAL*8 E2,EB,0B,EB2,EBS,PI,TW,A 


PI=3.141592654D+00 

TW=2.0D+00 

E2=EPS*EPS 

EB=1.0D+O0-E2 

OB=OMS-VPS*«TW 

EB2=EB*EB 

A=PAR/EB2 

EBS=EB**0.5 

P1=-Ax* (PI* (1.O0D+00+E2*TW) *VEP/EBS+OB*E2*TW) 
P2=Ax (EPS*VEP*TW*TW+OB*PI*EPS*EBS /TW) 
RETURN 

END 


Subroutine ROTTSA to calculate rotor angular 
travel, speed and acceleration 


SUBROUTINE ROTTSA(MAR,T,OMI,OMR,TA,OMT,OMS,OMA) 
REAL*8 T,OMI,OMR,TA,OMT,OMS,OMA,A,TX 
IF(MAR .EQ. 2) GO TO 2 

IF(MAR .EQ. 3) GO TO 3 

A=(OMR-OMI) /TA 

TX=T/TA 

OMA=A* (1. O0D+00-TX) *2.O0D+00 
OMS=A*Tx (2. OD+OO-TX) +OMI 

OMT=A*xT**2* (1.O0D+00-TX/3.0D+00) tOMI*T 
RETURN 

OMA=0.0D+00 

OMS=OMR 

OMT=OMR*T- (OMR-OMI) *TA/3.0D+00 

RETURN 

END 


Subroutine SPLINT to integrate discrete functions, 
using cubic spline method 


SUBROUTINE SPLINT(N,DX,X,Y,B,C,F1,F2,F3) 
REAL*8 DX,X(N),Y(N) ,B(N) ,C(N) ,F1,F2,F3,SUM1,SUM2,T 
NM1=N-1 
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K=1 

X(1)=0.0D+00 

DO 1 I=2,N 
X (1) =X (I-1) +DX 

C2 yay (2)-Y (1). /Dx 

DO 10 I=2,NM1 
B (I) =DX*4.0D+00 
e(i+1)= (Wy Gi +1)-Y G1)» /Dx 


CCT) =e (1¥1)\-c (1) 
B(1)=-DX 
B(N) =-DX 


C(1)=0.0D+00 
C(N) =0.0D+00 
IF (NiVEQ293) GO "TONS 
CC 1JsXCKS)RC(2))/ CDEX22 ODHOO) 
c(N) =(C (N- 1) -C:(N=-2) ) / (OX*2.0D+00) 
C(1)=C (1) *DX/3.0D+00 
C(N) =-C (N) *DX/3.0D+00 
DO 20 I=2,N 
T=DX/B C= 1) 
B(1I)=B (I) -T*DX 
e(1) =e (1) ete (i=1) 
C(N)=c (N) /B(N) 
DO 30 IB=1,NM1 
I=N-IB 
G(T He) zpxkxc (+1) )/B (1) 
SUM1=0.0D+00 
SUM2=0.0D+00 
DO 40 I=1,NM1 
SUM1=SUM1+Y (1) +Y (I+1) 
SUM2=SUM2+C (I) +C (I+1) 
T= (SUM 1-SUM2*DX**2/2.0OD+00) *DX/2.0D+00 
PECKeSEO. 1) F1= 3 
TECK NEO. 2) F2=T 
TECKICEQ. 3) F3=T 
DO 50 I=1,N 


¥ (i) =Y 1) *x (1) 
K=K+1 
IF(K .EQ. 4) GO TO 6 
GO TO 5 
RETURN 
END 
Subroutine to calculate components or Load Vector. 


Subroutine uses 4-point Gaussian quadrature to 
evaluate integrals. 


SUBROUTINE LOADYB(MAR,T,DT,TA,MR,OMR,OMI,ECC,TG,HG 
,F1,F2,F3) 
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REAL*8 T,DT,TA,MR,OMR,OMI,ECC,TG(4) ,HG(4) ,F1,F2,F3 
REAL*8 TX,OMT,OMS,OMA,FX,A,X,XT 


F1=0.0D+00 

F2=0.0D+00 

F3=0.0D+00 

IF(MAR .EQ. 2) GO TO 2 

IF(MAR .EQ. 3) GO TO 3 

A= (OMR-OMI) /TA 

DO 10 I=1,4 
TX=DT* (TG (I) +1.0D+00) /2.0D+00 
X=T+TX 
XT=X/TA 
OMT=Ax*xX*x*2%* (1. OD+00-XT/3.OD+00) tOMI*X 
OMS=A*X* (2. OD+00-XT) +OMI 
OMA=Ax (1.0D+O0-XT) *2.O0D+00 
FX=-MR*ECC* (OMS**2*DCOS (OMT) +OMA*DSIN (OMT) ) 
F1=F1+HG (1) *FX 
F2=F2+HG (1) *TX*FX 
F3=F3+HG (1) *«TX**2*FX 

GO TO 4 

A= (OMR-OMI) *TA/3.0OD+00 

DO 20° 1=154 
TX=DT* (TG (I) +1.0D+00) /2.0D+00 
X=T+TX 
OMT=OMR*X-A 
FX=—-MR*ECC*OMR**2*DCOS (OMT) 
F1=F1+HG (1) *FX 
F2=F2+HG (I) *TX*FX 
F3=F3+HG (I) *TX**2*FX 

CONTINUE 

A=DT/2.0D+00 

Fi=F1ix*xA 

F2=F2*A 

F3=F3*A 

RETURN 

END 
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